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Introduction 


Thrust  oscillations  generally  occur  in  large  solid  propellant  motor  during  a  few  seconds  of  the 
flight,  they  exactly  correspond  to  the  frequencies  of  the  cavity  modes.  The  excitation  of  cavity 
modes  has  been  reproduced  in  small  scale  solid  propellant  motors  and  also  in  cold  gas  facilities. 
However  as  extensively  exposed  by  F.  Vuillot  in  the  present  VKI  lecture  series,  there  are  different 
possible  causes  of  this  excitation  and  all  of  them  may  play  a  role.  Among  them,  as  suggested  by 
F.  Vuillot,  see  [1],  there  is  the  intrinsic  instability. 

Let  us  start  with  very  simple  models.  Let  be  C{uS)  a  linear  differential  operator  for  which 
some  values  of  to  correspond  to  eigenmodes.  It  means  that  for  these  values,  it  exists  a  non  zero 
solution  uu  (uu  ^  0)  such  as  :  £(u>)(uu)  =  0.  Now,  ui  is  assumed  to  be  such  an  eigenvalue  and 
Uuj  the  associated  eigenfunction.  If  the  system  is  forced  by  a  non  zero  right  hand  side  term  /, 
we  have  to  solve  :  C(uj)(v)  =  f.  Thanks  to  the  linear  nature  of  the  differential  operator,  if  / 
may  be  decomposed  in  different  terms  /  =  ^  /),  the  solution  v  is  the  exact  superposition  of  the 
solution  of  each  isolated  forcing  v  =  Y^,vi  with  C(u)(vi)  =  fi.  If  /,  does  not  correspond  to  the 
eigenmode  uu,  the  amplitude  Vi  remains  as  small  as  the  one  of  the  forcing  fi.  On  the  other  hand, 
if  fi, u>  is  an  eigenmode  (proportional  to  uw),  then  the  amplitude  of  Vi,u  may  become  very  large 
(it  is  the  so-called  secular  terms  in  the  framework  of  the  multiple  scale  analysis,  it  is  the  so-called 
resonance  phenomenon  generally  speaking).  Consequently,  after  some  transient,  the  contribution 
of  VitUJ  in  the  sum  v  =  Y^vi  will  be  dominant,  v  ~  Vi>u ,  :  the  system  has  selected  a  single  behaviour, 
the  one  corresponding  to  the  resonant  forcing. 

Coming  back  to  the  solid  propellant  motor,  the  previous  operator  £(w)  may  be  associated 
to  the  governing  problem  for  the  cavity  mode,  the  matter  is  to  determine  the  origin  of  the  forcing. 
An  intrinsic  instability  mechanism  has  been  suspected  by  Frangois  Vuillot  to  be  responsible  of 
the  forcing.  If  it  is  the  case,  the  forcing  would  be  itself  an  eigenmode,  but  for  another  dynamical 
system  ! 

The  present  course  is  limited  to  the  presentation  of  the  linear  stability  theory.  We  will  show 
the  general  philosophy  which  is  behind  a  so-called  “instability” ,  the  application  of  it  for  the  special 
flow  which  is  representative  of  the  one  in  a  solid  propellant  motor.  We  will  also  show  that  the 
present  theory  is  very  strange  from  a  theoretical  point  of  view,  but  very  efficient  in  comparison 
with  the  available  experimental  results.  Finally  we  will  conclude  as  a  “Public  Prosecutor”  against 
the  intrinsic  instability  modes  which  are  suspected  to  be  the  forcing  terms  with  respect  to  the 
emergence  of  the  dangerous  cavity  modes. 
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Chapter  1 

Geometry,  presentation 


1.1  Experimental  facilities 

In  real  motors,  the  conditions  are  very  severe  (two-phase  flow,  high  pressure,  high  temperature),  so 
that  detailed  measurements  inside  the  motor  are  nearly  impossible.  Cold  gas  experimental  facilities 
have  been  thus  especially  designed  for  fundamental  studies.  They  are  much  less  expensive  than 
the  experiment  which  use  solid  propellant  and  detailed  measurements  are  possible  ;  on  the  other 
hand  some  physical  effects  are  missing  such  as  the  reactive  two-phase  flow  or  the  slow  regression 
of  the  wall  (due  to  the  combustion  of  the  propellant) . 

For  cold  gas  experiment,  the  flow  inside  the  motor  which  comes  from  the  burning  surface  in 
real  motors  is  simulated  by  a  wall  injection  of  cold  gas  (air  for  example).  This  is  usually  realized 
by  using  a  porous  wall  constituted  by  small  bronze  particles  which  are  then  compacted. 

In  the  following,  two  types  of  geometry  will  be  considered. 

•  The  first  one  is  of  rectangular  type.  We  will  take  benefit  from  only  one  facility  of  this  kind  : 


Figure  1.1:  VECLA  facility 
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the  one  called  “VECLA”,  extensively  studied  by  G.  Avalon  in  ONERA  Palaiseau,  see  [2]  even 
if  other  facilities  such  as  the  so-called  “Micatl”  carried  out  in  Poitiers,  see  [3],  with  which  very 
interesting  results  have  been  obtained.  A  photo  of  VECLA  is  given  in  figure  1.1,  whereas  a 
sketch  is  given  in  figure  1.2.  Its  length  and  its  breadth  are  fixed,  603  mm  and  60  mm  respectively. 


Pc 

Tp 


Adjustable 

nozzle 


feeding  tube 


Figure  1.2:  Sketch  of  the  VECLA  facility 


On  the  other  hand,  the  height  can  be  modified,  10  mm,  20  mm  and  30  mm  have  been  tested.  Air 
is  injected  only  through  the  lower  wall.  The  upper  wall  is  solid  and  is  perforated  by  small  holes 
in  order  to  maintain  a  hot  wire  anemometer  used  for  measurement  of  the  instantaneous  velocity. 
The  two  lateral  walls  are  in  plexiglas  allowing  thus  flow  visualisations.  In  the  theory  presented 
in  the  following,  we  will  considerer  a  plane  case  for  which  both  lower  and  upper  walls  are  used 
for  the  air  injection.  Finally  the  system  is  closed  on  one  side  by  a  front  wall  and  ends  on  the 
other  side  by  either  nothing  or  by  a  non  injecting  nozzle  leading  to  a  throat.  In  the  operating 
conditions,  the  flow  is  sonic  at  the  throat.  Without  nozzle,  the  injection  velocity  may  be  varied 
continuously  during  the  experience.  It  can  also  be  fixed,  typical  values  are  around  1  m/s.  The 
nozzle  is  adjustable,  as  indicated  in  figure  1.2,  so  that  the  injection  velocity  may  be  changed  but 
not  continuously  during  an  experiment.  In  addition  to  some  technical  measurements  necessary 
for  the  control  of  the  flow,  the  measurements  consist  in  the  fluctuating  pressure  at  the  front 
wall  (denoted  by  Pc  in  figure  1.2)  and  the  velocity  (mean  value  and  fluctuating  one)  inside  the 
flow  by  using  a  hot  wire.  This  one  can  be  moved  at  several  distances  from  the  front  wall  and 
different  distances  from  the  porous  wall  can  be  analysed.  By  using  the  so-called  periodogram 
process,  the  temporal  velocity  signal  is  finally  converted  into  a  spectral  representation.  Very 
interesting  and  recent  results  obtained  with  VECLA  may  be  found  [4]. 

•  The  second  geometry  is  closer  than  the  one  of  real  solid  propellant  motors,  it  is  a  cylinder. 
Results  of  different  facilities  will  be  used  in  the  present  document.  The  oldest  one  (whose  results 
are  used  in  the  present  document)  is  located  in  the  United  States  and  has  been  extensively  tested 
by  Brown  and  co-workers,  see  [5]  for  example.  Another  one  is  in  Belgium  and  has  been  realized 
and  studied  by  Jerome  Anthoine,  see  [6],  manager  of  the  present  lecture  series  !  Finally  a  more 
recent  one  is  currently  investigated  at  ONERA-Palaiseau  by  G.  Avalon.  A  photo  of  this  last 
facility  is  provided  in  figure  1.3.  In  that  case,  there  is  a  unique  cylindrical  porous  wall,  through 
which  cold  gas  is  injected.  The  diameter  is  60  mm.  This  case  is  obviously  closer  to  the  real 
geometry,  but  only  one  diameter  can  be  used  for  a  given  system.  As  for  the  first  geometry, 
there  is  a  front  wall  and  either  a  free  exit  section  (as  it  is  the  case  for  the  photo  1.3)  or  a 
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Figure  1.3:  Photo  of  the  VALDO  facility 


throat.  As  with  VECLA,  for  a  free  exit  section,  the  injection  velocity  can  be  varied  more  or 
less  arbitrarily.  Front  wall  pressure  is  measured  and  a  hot  wire  is  also  used,  it  may  be  observed 
in  photo  1.3.  But  for  this  type  of  installation  with  a  hot  wire  passing  through  the  porous  wall, 
very  few  distances  from  the  front  wall  can  be  analysed.  In  order  to  obtain  finer  experimental 
results  for  the  evolution  of  the  velocity  with  respect  to  the  distance  from  the  front  wall,  the  hot 
wire  has  been  installed  on  a  very  fine  and  long  blade  passing  directly  through  the  exit  section. 


1.2  Notations 


1.2.1  Plane  case 

As  sketched  in  figure  1.4,  the  x  coordinate  is  defined  in  the  direction  perpendicular  to  the  front 


porous  wall 


streamline 


I— 4-^f— ■ 4—1 f— • -[—■ ]—■ 4— |- 


2h 


V, 


mj 


y‘ 

Z 

Figure  1.4:  Plane  notations 


wall,  x  =  0  being  located  at  this  wall.  The  coordinate  denoted  by  y  defines  the  distance  from  the 
upper  and  lower  walls,  y  =  0  is  located  in  the  symmetry  plane.  Finally,  the  coordinate  z  defines 
the  broadness.  The  distance  between  the  upper  and  lower  walls  is  denoted  by  2 h,  so  that  h  is  the 
physical  distance  between  the  two  walls  when  considering  the  VECLA  facility,  for  which  the  non 
porous  wall  is  located  at  y  =  0.  The  norm  of  the  injection  velocity  through  these  walls  is  assumed 
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to  be  constant  and  uniform,  it  is  noted  V)nj. 


1.2.2  Axisymmetric  case 

As  in  the  first  case,  x  represents  the  distance  from  the  front  wall.  Otherwise,  r  denotes  the  radial 
coordinate  and  9  the  azimuthal  one.  The  geometry  is  sketched  in  figure  1.5.  The  diameter  is  noted 
by  2 h  and  the  norm  of  the  injection  velocity  by  hjnj,  as  for  the  plane  case. 


porous  wall 


— >x  _ streamline _ 


2h 


V 


mj 


Figure  1.5:  Axisymmetric  notations 


1.2.3  Dimensionless  quantities,  Reynolds  number 

In  the  following,  all  the  quantities  are  dimensionless  according  to  the  reference  velocity  Idnj  and  the 
reference  length  h,  it  means  that  the  lengths  are  scaled  by  the  radius  (and  not  by  the  diameter)  of 
the  cylinder  in  case  of  the  axisymmetric  geometry.  The  reference  time  is  h/V\u j  and  the  pressure  is 
scaled  with  pI/A,  p  being  the  density  (assumed  to  be  constant).  Thus,  in  the  equations  it  appears 
a  Reynolds  number  based  on  the  injection  velocity  : 

v 

where  v  represents  the  kinematic  viscosity.  In  the  operating  conditions,  this  Reynolds  number  is 
of  order  1000.  In  some  papers,  like  [7],  another  Reynolds  number  is  used  which  is  based  on  the 
longitudinal  velocity  in  the  middle  of  the  channel  at  the  considered  distance  from  the  front  wall, 
see  also  figure  4.3  in  the  following. 

1.3  General  equations 

1.3.1  Used  assumptions 

As  previously  mentioned,  flow  injection  is  assumed  to  be  steady  and  uniform  (independent  of  the 
location)  and  strictly  perpendicular  to  the  porous  wall.  Concerning  the  theoretical  results,  the 
geometry  is  assumed  to  be  constant.  In  fact,  as  demonstrated  in  [8],  a  slow  regression  of  the 
porous  wall  (compatible  with  the  real  time  scale  of  the  motion  of  the  propellant  surface  in  real 
motors)  does  not  modify  significantly  the  results  presented  in  the  following. 

Finally  the  total  length  of  the  channel  is  assumed  to  be  sufficiently  short  so  that  the  flow 
remains  subsonic  inside  the  channel  (it  will  be  shown  below  that  the  flow  is  uniformly  accelerated 
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in  the  x  direction  due  to  the  injection).  As  the  temperature  field  is  numerically  observed  to  be 
more  or  less  constant,  we  also  assume  that  the  flow  is  incompressible.  Again  this  assumption  has 
been  validated,  see  [9]. 

1.3.2  Plane  geometry 

Assuming  that  the  flow  remains  two-dimensional,  all  the  quantities  depend  only  on  the  time  t  and 
on  the  spatial  coordinates  ( x ,  y).  The  two  components  of  the  instantaneous  flow  velocity  are  noted 
(u,v),  the  pressure  p,  so  that  the  Navier-Stokes  equations  written  for  the  incompressible  flow  in 
the  cartesian  (x,  y)  coordinate  system  are  : 


du 

dx 

du 

at 

dv 

at 


£>o 

dy 
_  du 


dx 

.  dv 
1  dx 


„du  dp  1 

Vdy  +  ai  =  Ra“ 


(1.1) 


.dv 

% 


f  =  Aas 

dy  1Z 


where  A  denotes  the  laplacian  operator.  The  boundary  equations  must  first  express  the  no-slip 
condition  together  with  the  injection  condition  at  the  two  horizontal  porous  walls.  Concerning  the 
front  wall,  the  no-slip  condition  leads  to  a  non  self-coherent  relation,  especially  at  the  two  corners 
(0,-1)  and  (0,1)  in  terms  of  (x,y)  coordinates.  Assuming  an  injection  at  these  points  is  not 
compatible  indeed  with  the  viscous  no-slip  relationship.  A  possible  solution  could  be  to  simulate 
a  boundary  layer  in  the  injection  process  in  the  neighbourhood  of  these  two  points.  A  simpler 
solution  however  consists  in  assuming  that  the  front  wall  acts  as  a  symmetry  plane  (allowing  a 
mirror  flow  for  negative  values  of  a;).  This  last  solution  leads  to  impose  that  only  the  u  component 
vanishes  at  x  =  0.  Finally  the  boundary  conditions  associated  to  the  equations  (1.1)  are  : 


Vye  [-1,+1]  u(0,y)  =  0 

\/x>0  u(x,  —  1)  =  0,  v(x,  —  1)  =  1,  u(x,  1)  =  0,  v(x,  1)  =  —  1 


(1.2) 


For  this  type  of  flow,  it  may  be  useful  to  work  with  a  stream  function  ip  related  to  the  velocity  by  : 


u  = 


dip 

dy 


v  =  — 


dip 

dx 


In  system  (1.1),  the  continuity  condition  is  automatically  satisfied  for  ip,  whereas  the  pressure  can 
be  eliminated  between  the  two  momentum  equations.  Differentiating  the  first  momentum  equation 
with  respect  to  y  and  the  second  one  with  respect  to  x  and  subtracting  then  these  two  equations 
lead  to  : 


d  f  d  _  d 
dt  \ayU~  dxV 


„  d  (  d  _  d  _ 

“s  y-s' 


_  a  t  s  _  o  \  i  . 

'’avU“-ai'V  =  Ka 


a  _  a  _ 

dyU  ~  dxV 


In  terms  of  stream  function,  this  writes  : 


d  ~  dip  d  -  dip  d  -  lAAr 
—  A^+  ——Ai’-  ——Aip=  —  AA  ip 
ot  oy  ox  ox  oy  7c 
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The  boundary  conditions  must  then  also  be  expressed  in  terms  of  ip  : 

dip 


Vy  e  [-1,+1] 


dy 


(0,y)  =  0 


Vx  >  0 


dip 

dy 


(*,-!)  =  0,  g(*,-l)  =  -l,  |(*,1)=0  g(x,l)  =  l 


(1.4) 


1.3.3  Cylindrical  geometry 

In  this  case,  p  remains  the  notation  for  the  pressure  and  the  three  components  of  the  instanta¬ 
neous  velocity  write  now  as  (ux,ur,ug),  in  the  cylindrical  ( x,r,0 )  coordinate  system  illustrated  in 
figure  1.5.  The  Navier-Stokes  equations  become  now  : 


dux  dur  ur 
dx  dr  r 

dux  _  dux 
^t+Ux^x+' 


dur 

~dt 


dug 

dt 


dur 


1  dug 

I - - 

r  dO 

dux 
r  dr 

dur 


=  0 

_  1  dux 


Ug 


r  dO 
1  dilr 


dp 

dx 


d2i 


d2i 


1Z  \  dx2 


dr 2 


1  dux 
r  dr 


1  d2ux 
r2  dO2 


u x  -X - r  ur  — —  +  Ug  -  — — 

dx  dr  r  dO 


d2i 


d2 


ur 


Ug_  +  dp 
r  dr 
1  dur  i 


~  dug  ,  _ 

^X  0  T  0 
dx  dr 


1Z  \  dx2  dr2  r  dr  r2 
dug  Ug  dug  UgUr  1  dp 


1  d2ur 
r2  dO2 


2  dug 
r2  dO 


(1.5) 


r  dO  r  r  dO 
1  f  d2ug  d2ue  1  due  _  Me 
7 Z  1  dx 2  dr2  r  dr  r2 


1  d2ug  2  dur 
r2  dO2  r2  dO 


Concerning  the  boundary  conditions,  the  same  previous  remark  about  the  injection  at  the  abscissa 
x  =  0,  close  to  the  front  wall,  can  be  expressed.  The  boundary  conditions  are  thus  : 


V<9<e[0,2tt[  Vr<E[0,l] 


t(O,r,0)=O 


Va;  >  0,  MO  e  [0,  2tt[  ux(x,1,0)=0,  ur(x,l,0)  = —1,  ug(x,l,0)  =  0 


(1.6) 


1.4  Basic  flow 


In  this  section,  we  will  first  determine  a  particular  steady  solution  of  the  Navier-Stokes  equations 
and  associated  boundary  conditions.  The  following  step  of  the  analysis  will  be  to  determine  its 
stability,  this  will  be  done  in  the  following  chapters.  The  physical  quantities  associated  to  the 
steady  solution,  we  are  looking  for,  will  be  noted  in  capital  and  overlined  letters. 


1.4.1  Plane  case 

It  is  possible  to  find  a  mathematical  steady  solution  of  the  equations.  A  self-similar  solution  firstly 
proposed  by  Berman,  see  [10],  may  be  sought  for  the  stream  function  of  the  mean  flow  : 

i>  =  xF(y) 
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that  is  U  =  xF'{y )  and  V 


—F  in  terms  of  the  velocity  components.  Equation  (1.3)  leads  to 


p'  p"  —  FF'" 


LpiIV) 

TZ 


associated  to  the  boundary  conditions  : 


(1.7) 


F(— 1)  =  -1  F'(-1)=0  F(l)  =  1  F'(1)=0 


(the  boundary  condition  at  the  front  wall  is  automatically  satisfied  with  this  self-similar  form). 
For  large  Reynolds  numbers,  a  very  accurate  approximation  of  the  basic  flow  is  obtained  with  the 
so-called  Taylor  solution,  see  [11]  : 


-  7tx  Try 
U  =  —  cos  — 
2  2 


f>  •  ny 

V  =  —  sin  — 
2 


(1.8) 


which  is  strictly  valid  for  the  equations  written  for  an  inviscid  flow.  The  (inviscid)  pressure  is  then 
given  by  : 

1 


P=  - 


(.  nyy 

V  2  J 


Pn 


where  Pq  is  a  constant.  The  associated  stream  function  is  : 


with  k  a  constant. 


ip  =  x  sin 


Try 

2 


+  k 


1.4.2  Axisymmetric  case 


In  the  cylindrical  case,  a  steady  axisymmetric  solution  may  be  also  determined  in  the  self-similar 
form  given  by  : 

1  (94*  -  1 94* 

Ux  =  Ur  =  — — —  Ug  =  0  with  4>(ar,r)  =  xf{r) 

r  or  r  ox 

where  4>  is  a  cylindrical  stream  function  of  the  flow.  The  function  /  satisfies  the  following  differ¬ 
ential  problem  with  the  prime  corresponding  to  the  differentiation  with  respect  to  r  : 


1 

77 


-  <  r  — 


r 


(0)  =  0 


-(0)=0  /'(1)=0 
r 


/ 


=  0 


/( 1)  =  0 


The  two  first  boundary  equations  written  above  impose  the  leading  order  of  the  behaviour  of  the 
function  /  close  to  the  axis  r  =  0  in  case  of  a  viscous  flow.  Assuming  a  regular  Taylor  expansion 
of  the  form  f(r)  =  ao  +  a\r  +  a,2r2 /2  +  . . .,  it  can  be  readily  proved  that  /  satisfies  f(r)  =  0(r3) 
for  r  close  to  0.  Without  viscosity,  only  the  second  boundary  condition  in  r  =  0  is  necessary  and 
the  leading  order  becomes  then  /(r)  =  0(7'2)  for  r  close  to  0. 

Again  it  is  possible  to  find  an  analytical  solution,  which  exactly  satisfies  the  inviscid  equa¬ 
tions  : 

9  -i  9 

_  ttt  -  1  nr  - 

Ux  =  7tx  cos  ^  Ur  =  —  sin  ^  Ug  =  0  (1.9) 
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with  the  pressure  given  by  : 


P=- 


7 T2X2 


2 


1 

^2 


+  Po 


where  Pq  is  a  constant.  This  solution  constitutes  a  good  approximation  of  the  viscous  case  for 
large  Reynolds  numbers  (typically  TZ  >  1000). 


1.4.3  Some  remarks 

Before  analysing  the  stability  properties  of  the  mean  flow  described  above,  it  may  be  interesting 
to  note  some  remarkable  features. 

•  First  the  approximate  analytical  solution  has  been  checked  to  be  quasi  superposed  to  the  exact 
self  similar  solution  for  the  two  types  of  geometry  at  least  for  the  Reynolds  numbers  range 
considered  in  the  present  document,  see  [8]  and  [12]. 

•  The  mean  flow  depends  on  two  spatial  coordinates  and  accordingly  there  is  a  non  zero  vertical 
(or  radial)  velocity  component,  so  that  the  flow  is  said  to  be  non  parallel.  The  plane  Poiseuille 
flow,  solution  in  an  infinitely  long  and  broad  non  injected  rectangular  duct,  is  strictly  parallel, 
there  is  only  one  spatial  direction  in  which  the  velocity  is  not  homogeneous  (the  distance  to 
the  wall)  and  consequently  the  streamlines  are  strictly  parallel.  Conversely,  the  Taylor  flow  is 
non  parallel,  the  streamlines  start  at  the  porous  lines,  turn  and  become  more  or  less  parallel  far 
from  the  front  wall.  This  is  illustrated  in  figure  1.6  for  both  types  of  geometry.  In  fact,  even  for 


Figure  1.6:  (Computed)  streamlines  of  the  Taylor  flow  (plane  geometry  in  the  upper  figure,  ax- 
isymmetric  geometry  in  the  lower  one). 

large  values  of  x,  it  remains  obviously  a  small  region  close  to  the  porous  wall  in  which  the  flow 
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is  highly  non  parallel.  This  feature  induces  some  difficulties  with  respect  to  the  mathematical 
form  of  the  perturbation  considered  in  the  stability  analysis.  This  point  will  be  discussed  in 
some  details  in  the  following,  see  section  4.1. 

•  This  x  dependence  describes  explicitly  a  linear  increase  of  the  longitudinal  component  of  the 
velocity. 

•  Comparing  the  two  types  of  geometry  in  figure  1.6,  it  appears  clearly  that  at  a  given  x  position 
the  curvature  of  the  streamlines  is  much  larger  in  the  axisymmetric  case  than  in  the  plane  case. 
This  will  strongly  affect  the  stability  results. 
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Chapter  2 

Linear  Stability  Theory 


2.1  A  short  philosophical  escape 

It  must  be  emphasized  that  the  mean  flow  described  in  the  previous  chapter  is  a  possible  solution  of 
the  equations.  However  finding  such  a  possible  solution  does  not  prove  that  this  solution  will  be  the 
observed  one  in  practice.  This  last  point  depends  in  fact  on  the  stability  of  the  proposed  solution. 
A  similar  problem  arise  for  example  with  a  small  ball  placed  on  the  top  of  a  large  sphere  (where 
“top”  refers  to  the  direction  of  the  gravity).  The  position  on  the  top  is  a  possible  equilibrium 
location,  however  as  it  is  unstable,  it  cannot  be  practically  observed. 

Coming  back  to  the  Taylor  flow,  the  same  question  is  relevant  :  is  the  particular  flow 
described  by  expression  (1.8)  or  (1.9)  observed  in  practice,  that  is  for  example  in  the  channel 
VECLA  ?  In  fact  the  answer  depends  on  the  injection  conditions  (injection  velocity  and  porosity 
of  the  material) .  The  injected  flow  probably  contains  some  turbulent  structures  which  may  excite 
continuously  the  channel  flow.  It  has  been  observed  that  for  a  given  porosity  if  the  injection  velocity 
is  too  large,  the  Taylor  flow,  as  it  is  described  by  expression  (1.8)  or  (1.9)  is  never  observed.  In  that 
conditions,  the  flow  inside  the  channel  is  turbulent  everywhere,  so  that  it  exists  large  turbulent 
fluctuations,  which  interact  between  themselves  leading  by  nonlinear  mechanisms  to  a  modification 
of  the  mean  flow.  The  Taylor  flow  which  is  a  possible  solution  of  the  averaged  equations  does  not 
correspond  then  to  the  averaged  observed  (or  measured)  flow.  This  case  will  be  no  longer  examined 
in  the  following  part  of  the  present  document. 

For  small  injection  velocities  and  for  small  values  of  the  porosity,  the  Taylor  flow  can  be 
observed,  at  least  in  the  upstream  part  of  the  channel.  In  fact,  the  measured  flow  does  not 
correspond  exactly  to  the  Taylor  flow,  it  always  exists  some  fluctuations  (for  example  those  induced 
directly  from  the  non-perfect  injection  system)  superimposed  to  this  theoretical  solution.  As  all 
dynamical  systems,  the  main  flow  which  is  continuously  excited  by  the  injection  system  exhibits 
two  types  of  response  with  respect  to  this  forcing.  First  there  is  the  so-called  forced  response,  whose 
amplitude  is  of  the  same  order  of  the  forcing  amplitude.  Secondly  it  may  exist  an  eigenresponse 
the  amplitude  of  which  can  be  dangerously  larger  than  the  one  of  the  forcing. 

The  amplitude  of  the  eigenresponse  may  be  in  fact  extremely  large,  it  means  that  a  micro 
phenomenon  (small  inhomogeneity  of  the  injection)  can  generate  a  macro  feature.  This  particular 
behaviour  is  the  signature  of  an  intrinsic  instability.  Obviously  determining  the  physical  character- 
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istics  of  the  eigenresponse  is  very  important  (amplification  rate,  frequency,  spatial  shape),  this  is 
the  goal  of  the  stability  studies.  A  well-known  and  classical  presentation  of  the  stability  theory  is 
published  in  [13],  whereas  a  general  presentation  of  the  practical  use  of  it  in  case  of  the  boundary 
layer  which  exhibits  some  common  feature  with  the  present  basic  flow  is  nicely  presented  in  [14] 
and  finally  the  first  attempt  of  a  theoretical  application  for  the  Taylor  flow  has  been  done  in  Russia 
by  Varapaev  and  Yagodkin  and  has  been  translated  in  [7]. 

2.2  Small  perturbation  technique 

According  to  the  previous  comments,  a  particular  solution  (the  Taylor  flow  in  the  present  case) 
is  supposed  to  be  known,  its  stability  is  going  to  be  analysed.  The  instantaneous  flow  is  then 
assumed  to  result  from  a  pure  superposition  of  this  particular  flow,  which  will  be  called  the  “basic 
flow”  in  the  following,  and  of  a  fluctuation  to  be  determined.  This  is  mathematically  imposed  by 
writing  : 

q  =  Q  +  q  (2.1) 

for  any  physical  quantity  q  (components  of  the  velocity,  pressure  and  also  for  example  temperature, 
mass  flux  density  in  case  of  compressible  fluids).  The  instantaneous  flow  is  assumed  to  be  realistic, 
thus  it  satisfies  the  general  governing  equations.  The  key  point  for  the  stability  analysis  is  the 
following  one.  As  we  only  focus  on  an  eigenresponse  (and  not  on  a  forced  response),  the  boundary 
conditions  written  for  the  instantaneous  flow  must  be  strictly  identical  to  those  imposed  for  the 
determination  of  the  basic  flow. 

Before  writing  the  complete  equations,  it  may  be  instructive  to  describe  formally  the  pro¬ 
cedure.  Let  us  represent  the  complete  equations  together  with  the  boundary  conditions  by  a  (non 
linear)  operator  £,  for  example  the  one  corresponding  to  the  Navier-Stokes  equations  and  associ¬ 
ated  boundary  conditions.  As  explained  above  both  basic  flow  and  instantaneous  flow  are  assumed 
to  satisfy  these  equations  and  boundary  conditions  : 

C{q)  =  0  and  C(Q)  =  0  (2.2) 

According  to  the  superposition  (2.1),  the  first  equation  becomes  :  C(Q  +  q)  =0.  Then,  the 
operator  C  being  regularly  dependent  on  the  physical  quantities,  a  first  order  Taylor  expansion  can 
be  written  : 

C(Q  +  q)  =  C(Q)  +M(Q).q  +  NLT 

where  NLT  means  nonlinear  terms  with  respect  to  the  fluctuating  quantities  q  and  A i(Q)  is  a 
linear  operator  function  of  the  basic  flow  which  applies  to  the  fluctuating  quantities.  It  represents 
formally  the  gradient  operator  of  C.  In  this  course,  the  fluctuation  is  assumed  to  be  small  in 
comparison  to  the  basic  flow  :  g<Q,  so  that  the  nonlinear  terms  are  small  in  comparison  to  the 
linear  ones.  Indeed  in  the  linear  stability  theory  these  nonlinear  terms  are  simply  dropped,  the 
previous  Taylor  expansion  is  thus  simplified  into  : 

C(Q  +  q)  ~  C{Q)  +  M(Q).q 

Finally,  as  both  instantaneous  flow  and  basic  flow  satisfy  the  equations  as  expressed  by  (2.2), 
previous  equation  becomes  : 

M(Q).q  =  0  (2.3) 
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which  represents  the  linearised  Navier-Stokes  equations  (and  boundary  conditions).  In  the  general 
case,  the  operator  A4(Q)  is  invertible  so  that  the  only  solution  is  q  =  0,  it  means  that  there  is  no 
eigenresponse.  However  is  some  very  specific  cases,  as  for  the  present  Taylor  flow,  possible  non 
zero  eigensolution  may  be  found.  This  is  mainly  the  case  when  the  basic  flow  exhibits  symmetry 
properties  (usually  invariance  with  respect  to  spatial  and  temporal  coordinates)  inducing  then  a 
special  mathematical  form  for  the  fluctuation.  This  is  described  in  the  next  section. 

2.3  Normal  mode  form 

In  this  section,  the  basic  flow  is  assumed  to  be  dependent  on  only  one  spatial  variable  :  y  in 
the  plane  case  and  r  in  the  axisymmetric  case.  This  is  the  parallel  assumption.  Of  course,  as 
mentioned  above,  this  is  not  strictly  the  case  for  the  Taylor  flow,  see  section  4.1.  In  fact,  to  be 
more  precise,  the  parallel  assumption  consists  in  assuming  that  only  the  physical  quantities  Q, 
which  are  written  in  the  linearised  operator  (2.3),  are  dependent  on  one  spatial  variable.  For 
example  for  the  plane  Poiseuille  flow,  which  is  strictly  parallel,  the  pressure  depends  on  x  too,  but, 
at  least  in  the  incompressible  approximation,  this  pressure  term  does  not  appear  explicitly  in  the 
linearised  equations  (2.3). 

Just  to  fix  the  ideas,  we  will  note  y  the  spatial  variable  in  the  non  homogeneous  direction  : 
Q  =  Q{y)-  Then,  equation  (2.3)  represents  a  system  of  partial  differential  equations  whose  coef¬ 
ficients  depend  on  the  basic  flow  quantities  Q ,  hence  on  y  only.  Thus  with  respect  to  each  of  the 
other  variables,  equation  (2.3)  appears  as  a  system  of  linear  equations  with  constant  coefficients. 
It  is  well  known  that  solutions  in  that  case  may  be  sought  in  an  exponential  form.  Following 
this  mathematical  result,  all  the  fluctuating  quantities  are  sought  with  the  so-called  normal  mode 
form  : 

q(x,  y,  t)  =  q(y)e^ax  ~  plane  case 

(2-4) 

q( x,  r,  9 ,  t)  =  q(r)e^ax  +  axisymmetric  case 

This  notation  needs  to  be  explained.  First,  some  quantities  have  complex  values,  for  example  i 
is  the  imaginary  number  satisfying  i2  =  —1,  whereas  the  physical  fluctuation  must  be  obviously 
real.  In  fact,  the  physical  fluctuating  quantities  correspond  to  the  real  part  of  the  right  hand  side 
terms  in  (2.4).  The  functions  q(y)  and  q(r)  are  complex  and  are  called  the  amplitude  functions, 
m  is  an  integer  and  is  the  azimuthal  wave  number  and  a  and  w  are  generally  complex  numbers. 
Introducing  real  part  and  imaginary  part  of  these  two  numbers  by  subscripts  “r”  and  “i”,  the 
fluctuating  quantities  q  write  : 

q( x,  y,  t)  =  q{y)e~aiX  +  Uit  el(arX  ~  ^ 

The  second  exponential  term  is  of  norm  1,  it  describes  thus  the  wavy  nature  of  the  solution  for  the 
fluctuation,  ar  is  the  longitudinal  wave  number,  ur  the  circular  frequency  with  /  =  ui/ 2n  being 
the  frequency  itself.  The  first  exponential  term  is  real,  it  describes  a  possible  amplification  of  the 
fluctuation  with  respect  to  the  time  and/or  with  the  distance  x  according  to  the  sign  of  u>i  and  a*. 
Usually  at  this  step,  two  types  of  theory  are  distinguished  :  the  temporal  theory  for  which  ati  =  0, 
the  fluctuations  only  grow  with  time  and  the  growth  is  governed  by  the  temporal  growth  rate  u>i 
and  the  spatial  theory  for  which  on  the  other  hand  we  have  uii  =  0,  the  fluctuations  only  grow  in  x 
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and  the  growth  is  governed  by  the  spatial  growth  rate  —a,.  Between  these  two  theories,  it  may 
be  interesting  to  recall  that  in  fact  both  theories  usually  provide  very  similar  results  by  using  the 
so-called  Gaster  transformation,  see  [15]  : 


—  with 

V0 


V  =  — 
9  da 


where  Vg  is  the  group  velocity  (such  an  approximation  is  strictly  valid  in  the  vicinity  of  the  neutral 
curve) . 

As  suggested  by  the  experimental  results  for  the  Taylor  flow,  which  will  be  given  in  sec¬ 
tion  3.4,  the  relevant  approach  in  the  present  case  is  the  spatial  theory.  The  mathematical  form 
of  the  perturbation  used  in  the  present  document  is  thus  : 


q(x,y,t)  =  q(y)e  aix  el(arx  urt)  _  q(y)el{ax  cot)  (2.5) 


with  ui  a  real  number  and  a  a  complex  one.  Assuming  that  the  perturbation  is  convected  in  the 
same  direction  as  the  basic  flow,  that  is  for  increasing  values  of  x,  the  spatial  evolution  of  the 
fluctuation  and  consequently  the  stability  of  the  basic  flow  depends  on  the  sign  of  cti,  the  results 
are  given  in  table  2.1.  The  same  conclusions  are  true  for  the  perturbation  form  in  the  axisymmetric 


sign  of  ccj  basic  flow 

ai  >  0  stable 

at  <  0  unstable 

Oj  =  0  neutral  (or  marginal) 

Table  2.1:  Stability  criterion  in  term  of  spatial  amplification  growth  rate 


geometry. 

Finally,  the  perturbation  given  by  relation  (2.4)  in  the  “plane  case”  corresponds  to  a  plane 
perturbation.  A  more  general  form  of  the  perturbation  must  include  the  z-dependence,  the  per¬ 
turbation  writing  then  : 

q(x,y,z,t)  =  q(y)ei(ax  +  Pz-“t) 

It  is  known,  see  [13]  for  example,  that  for  a  strictly  parallel  basic  flow,  it  is  sufficient  to  limit 
the  stability  analysis  to  two-dimensional  modes,  as  written  in  (2.4).  This  result  comes  from  the 
so-called  Squire’s  theorem.  This  one  cannot  be  theoretically  applied  due  to  the  nonparallel  terms 
which  are  kept  in  the  stability  equations.  However,  systematic  stability  computations  performed 
for  a  three-dimensional  perturbation  showed  that  the  most  amplified  mode  (which  corresponds  to 
the  largest  growth  rate  — a i)  is  obtained  for  (3  =  0,  that  is  for  a  two-dimensional  mode.  This 
explains  why  the  perturbation  is  assumed  to  be  in  the  form  (2.4)  in  the  following. 

2.4  Dispersion  relation 

The  problem  now  is  to  determine  the  amplitude  functions,  the  frequency  and  of  course  the  complex 
wave  number  a.  The  perturbation  is  required  to  satisfy  the  linearized  equations  (2.3).  Conse¬ 
quently,  the  form  (2.5)  is  introduced  into  these  linearized  equations,  this  yields  : 

A f(Q,  a,u).q  =  0 
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where  A f  is  now  an  ordinary  differential  operator  with  respect  to  y  only.  Differentiations  with 
respect  to  t  and  to  x  simply  become  multiplications  by  —iui  and  by  ia  respectively.  In  order  to 
find  a  non  zero  solution  q ,  it  is  necessary  to  find  values  of  (a,  u>)  such  as  A f(Q,  a,  u>)  is  non  invertible. 
The  mean  flow  depends  on  x  and  on  y  but  practically  the  stability  analysis  is  performed  for  a  fixed 
value  of  x.  On  the  other  hand,  as  this  operator  corresponds  to  the  linearized  form  of  the  Navier- 
Stokes  equations  (1.1),  the  Reynolds  number  TZ  is  a  parameter.  Finally,  the  formal  linearized 
equations  may  be  written  as  : 

Af(x,  TZ,  a,  ui).q  =  0  (2-6) 

Now  we  have  to  find  conditions  for  which  J\f(x,TZ,a,uj)  is  not  invertible.  Again,  it  can  be  first 
noted  that  this  operator  is  usually  invertible,  however  it  may  be  possible  that  for  some  specific 
conditions  between  the  parameters  (x,7 Z,a,co),  the  operator  is  not  invertible.  These  (not  known 
and  not  obvious)  specific  conditions  are  noted  by  a  relation  : 

ZF(x,TZ,a,cu)  =  0  (2-7) 

which  is  called  a  dispersion  relation,  because  for  x  and  7 Z  being  fixed,  this  equation  binds  the 
frequency  and  the  phase  velocity  u)/a.  Of  course,  this  equation  is  usually  not  explicit,  numerical 
computations  are  necessary.  They  are  often  performed  in  two  steps  (eigenvalues  and  then  eigen¬ 
functions).  First,  the  differential  continuous  problem  being  discretised,  a  complex  eigenvalue  (a 
for  example)  is  searched  in  such  a  way  that  a  non  zero  solution  for  the  perturbation  may  exist. 
In  this  step,  equation  (2.7)  is  solved  in  fact,  even  if  the  latter  is  not  written  explicitly.  Then  the 
optional  second  step  consists  in  determining  the  non  zero  solution,  with  the  obtained  eigenvalue 
(which  allows  a  non  zero  solution).  Namely,  this  non  zero  solution  is  the  eigenfunction  associated 
to  the  obtained  eigenvalue.  Some  insight  about  the  numerical  aspects  are  given  in  appendix. 


2.5  Linearised  equations 

2.5.1  Plane  case 


As  explained  before,  the  small  perturbation  technique  is  used,  by  decomposing  each  quantity 
as  written  in  relation  (2.1).  In  this  section,  the  perturbation  is  written  with  the  normal  mode 
form  (2.4)  even  if  it  is  not  justified,  according  to  the  x  dependence  of  the  mean  flow.  This  problem 
related  to  the  non  parallel  effects  will  be  considered  in  section  4.1. 

For  the  plane  case,  it  is  possible  to  use  either  the  primitive  variables  ( u,v,p )  or  the  stream 
function.  Let  us  start  with  the  first  formulation.  The  procedure  given  in  section  2.2  leads  with  (1.1) 
to  : 

.  ,  dv 
lau  +  — —  =  0 
dy 


—iuu  +  iaUu  +  u— - 1-  V  — — I-  v— - K  iup  =  —  [  —  —  a2u 


— iuv  +  iaUv  + 


.dU 

dx 

dV 
L  dx 


^du 

dy 

dv 


V 


.dU 

dy 

dV 


1 

n 


dy  dy 


dp 

dy 


d2u 
dy2 

1  (d2v  2, 
TZ  [jdy2  ~  a  1 


This  system  is  associated  to  the  following  boundary  conditions,  coming  from  (1.2)  : 


(2.8) 


«(— 1)  =  v(— 1)  =  u(l)  =  D(l)  =  0 
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It  can  be  easily  shown  that  system  (2.8)  may  be  written  under  the  form  of  a  first  order  ordinary 
differential  linear  problem  applied  to  the  vector  (ft,  du/dy ,  v,p)  of  dimension  4,  so  that  this  system, 
associated  to  the  four  boundary  conditions  above  is  a  well-posed  problem.  It  remains  to  solve  it 
numerically,  see  the  appendix  for  some  numerical  details. 

The  stream  function  formulation  has  the  advantage  to  define  only  one  scalar  unknown  (the 
stream  function  amplitude).  The  superposition  (2.1)  is  introduced  in  equation  (1.3).  After  lineari¬ 
sation,  this  leads  to 

/  1  \  ,  /)2t>  ^  pp.rj  n 

-iu>  +  iaU  +  VD  —  — (-D2  -  a2)  ( D 2  -  a2)tp  -  -  ia -—.ip  =  0  (2.9) 

\  K  )  dy2  dyA 

where  D  represents  the  differentiation  with  respect  to  y.  This  equation  corresponds  exactly  to 
equation  (2.3)  in  case  of  equation  (1.3).  Associated  to  this  equation,  the  following  boundary 
conditions 

^(-1)  =  Dtp(- 1)  =  ^(1)  =  Dtp(  1)  =  0  (2.10) 

coming  from  (1.4),  must  be  satisfied.  The  stream  function  of  the  perturbation  is  thus  required 
to  satisfy  again  a  fourth  order  homogeneous  ordinary  differential  equation  (2.9)  associated  to  four 
homogeneous  boundary  conditions.  Obviously,  the  trivial  if)  =  0  remains  a  possible  solution.  The 
goal  is  to  find  if  for  specific  conditions  between  a  and  lo  (dispersion  relation),  another  solution 
may  exist.  It  is  also  important  to  note  that  equation  (2.9)  contains  three  supplementary  terms 
(due  to  the  non  zero  V  factor)  in  comparison  with  the  classical  Orr-Sommerfelcl  equation.  The 
latter  is  obtained  as  explained  before  but  for  a  strictly  parallel  mean  flow,  that  is  with  a  velocity 
of  the  form  (U(y),0).  These  additional  terms  are  those  associated  with  an  odd  order  of  derivation 
in  (2.9). 


2.5.2  Axisymmetric  case 


The  same  procedure  can  be  applied  in  the  axisymmetric  geometry.  Like  the  plane  case,  the  basic 
flow  does  not  depend  on  r  only,  so  that  the  normal  mode  is  theoretically  not  applicable.  However 
using  it  leads  to  : 


.  „  dur  ur  im  „ 

iaux  +  — - 1 - 1 - ug  =  0 

or  r  r 


.  „  .  -  dUx  „  -  dux  dUx  „  .  „ 

-iuux  +  iaUxux  +  -w—Ux  +  Ur— - 1-  -7—ur  +  lap 

ox  or  or 


1 

=  I  —a~ur 

n 


d2ux  1  dur 


dr 2  r  dr  r2  3 


.  „  -  dUr  „  -  dur  dur  „  dp 

- lUJUr  +  iaUxUr  +  — — Ux  +  t/r— - 1-  — —  Ur  +  — 

dx  dr  dr  dr 


1 

=  I  —a~ur 

n 


d2ur  1  dur 
dr2  ^  r  dr 


m 


2 im  , 


9  9  ur 


-Ug 


■  -  ,  •  fr  -  ,  fr  9ue  .  Uriig  im  „ 

—iumg  +  iaUxug  +  Ur— - 1 - 1 - p 

dr  r  r 


d2iig  1  dug  Ug  TO2 


2 im  , 


dr2 


dr 


- v - TrUg 


(2.11) 
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This  system  of  linear  ordinary  differential  equations  is  associated  to  the  boundary  conditions  : 

«x(l)  =  Ur{  1)  =  'fifl(l)  =  0 

expressing,  as  for  the  plane  case,  that  there  is  no  fluctuating  velocity  (for  determining  the  eigen- 
response)  at  the  porous  wall.  The  numerical  resolution  of  system  (2.11)  presents  a  (small)  difficulty 
at  the  axis  r  =  0,  as  usual  when  using  the  cylindrical  coordinates  system,  see  [16]  for  a  possible 
treatment  of  this  singularity. 


RTO-EN-023 


8-23 


Motor  Flow  Instabilities  -  Part  2: 

Intrinsic  Linear  Stability  of  the  Flow  Induced  by  Wall  Injection 


ORGANIZATION 


8-24 


RTO-EN-023 


Motor  Flow  Instabilities  -  Part  2: 
Intrinsic  Linear  Stability  of  the  Flow  Induced  by  Wall  Injection 


Chapter  3 

Stability  results 


3.1  Eigenmodes 

The  first  task  is  to  determine  the  solutions  of  the  dispersion  relation  (2.7)  if  any.  As  the  conclu¬ 
sions  are  very  similar  in  the  other  cases,  we  will  only  consider  the  plane  case  formulated  by  the 
stream  function,  see  (2.9).  Then,  in  order  to  obtain  the  “complete”  spectrum  (i.e.  the  set  of  the 
eigenvalues)  of  the  linearised  equations,  it  is  more  convenient  to  fix  a  and  to  search  all  possible 
solutions  in  terms  of  u>  for  TZ  and  x  fixed.  The  reason  comes  from  the  fact  that  u>  appears  linearly 
in  the  Navier-Stokes  equations,  and  a  non  linearly  because  it  exists  terms  in  d2 /dx2  which  lead  to 
the  factor  a2.  More  precisely,  equation  (2.9)  may  be  rewritten  in  : 

|  (iaU  +  VD-  i( D 2  -  a2))  (D2  -  a2)  -  -  ia^  J  */>  =  itu(D2  -  a 2)^ 

which  is  formally  of  the  type  Ai[>  =  uiBip,  that  is  of  the  type  of  a  standard  (generalised)  eigen¬ 
value  problem.  The  differential  problem  being  discretised,  many  mathematical  libraries  allow  the 
determination  of  the  “complete”  spectrum. 

As  an  example,  the  following  parameters  have  been  fixed  :  TZ  =  1000,  x  =  10  and  a  =  4. 
Some  eigenvalues  are  plotted  in  figure  3.1  in  the  complex  w-plane.  The  horizontal  line  u>i  =  0 
separates  the  instability  zone  from  the  stability  zone.  The  higher  the  points  are,  the  most  dangerous 
they  are  (according  to  the  temporal  theory).  In  the  instability  region,  only  two  amplified  modes 
have  been  found  :  ui  £  {31.59  +  1.6946*,  31.622  +  1.6296*}.  They  are  thus  very  close  together  and 
seem  in  the  figure  to  be  only  one.  The  other  modes  are  damped  modes,  at  least  for  the  chosen 
values  of  the  parameters.  Moreover  for  smaller  value  of  u it  exists  many  other  modes,  including 
possible  continuous  branches. 

To  analyse  the  two  amplified  modes  and  their  differences,  it  is  generally  fruitful  to  calculate 
the  corresponding  eigenfunctions.  Moreover,  in  order  to  illustrate  the  strong  correspondence  be¬ 
tween  the  temporal  and  the  spatial  theory,  we  now  move  to  the  spatial  theory  and  fix  cu  to  a  real 
number  by  u>  =  31.6,  in  the  range  of  the  two  previous  modes.  Then  by  using  a  shooting  method 
(Newton  convergence  from  an  initial  guess  of  the  root),  two  amplified  modes  are  obtained,  as  in 
the  temporal  theory.  They  are  given  by  a  =  4.0037  —  *0.366  and  a  =  4.016  —  *0.373.  The  group 
velocity  may  be  calculated,  the  nearly  same  value  is  found  for  both  modes  Vg  —  4.4.  Gaster’s 
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Figure  3.1:  Spectrum,  1Z  =  1000,  x  =  10,  a  =  4,  300  points  in  the  collocation  method 


transformation  provides  then  temporal  amplification  rates  which  are  of  order  1.66,  they  constitute 
thus  good  predictions  of  the  temporal  growth  rates  which  have  been  computed  directly  in  the 
temporal  theory. 

The  corresponding  eigenfunctions  are  plotted  in  figure  3.2  for  the  first  one  and  in  figure  3.3 
for  the  second  one.  The  main  difference  between  these  two  modes  occurs  close  to  the  axis  y  =  0. 
For  the  first  mode,  plotted  in  figure  3.2,  the  transverse  velocity  v  is  zero  on  the  axis,  whereas  it 
is  the  longitudinal  velocity  u  which  vanishes  at  the  axis  for  the  second  mode.  The  first  mode  is 
called  a  varicose  mode,  and  the  second  one  a  sinuous  mode.  Concerning  the  facility  VECLA,  it 
can  be  noted  that  with  the  solid  wall  (non  porous)  placed  in  the  upper  limit  on  the  duct,  only 
the  varicose  modes  can  be  observed,  the  sinuous  ones  are  not  compatible  with  a  non-penetration 
condition  at  this  upper  wall. 


3.2  Amplitude  and  n  factor 

3.2.1  Definition  of  the  n  factor 

In  agreement  with  the  experimental  results,  see  section  3.4  for  validations  and  explanations,  the 
relevant  theory  is  the  spatial  one,  so  that  u>  is  a  real  number  characterising  the  frequency  of  the 
wave.  The  dispersion  relation  is  thus  solved  for  fixed  values  of  (77,  x,  u>)  (and  additionally  with  the 
azimuthal  wave  number  m  being  fixed  for  the  axisymmetric  geometry)  and  the  complex  a  value  is 
computed  in  order  to  satisfy  the  dispersion  relation  (2.7). 

The  opposite  of  the  imaginary  part  an  of  the  complex  wave  number  is  the  growth  rate  in 
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Figure  3.2:  Norm  of  the  longitudinal  and  transverse  velocity  component  of  the  eigenmode  a  = 
4.0037  —  *0.366  obtained  for  1Z  =  900,  x  =  10  and  u>  =  31.6 


Figure  3.3:  Norm  of  the  longitudinal  and  transverse  velocity  component  of  the  eigenmode  a  = 
4.016  —  z0.373  obtained  for  7 Z  =  900,  x  =  10  and  w  =  31.6 
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the  x  direction.  As  it  has  been  developed  for  the  boundary  layer  [17],  it  is  thus  possible  to  define 
an  amplitude  A  for  the  eigenmode  by  : 


1  dA 
~ai  ~  A~fa 


As  suggested  by  (2.7),  the  growth  rate  depends  on  the  parameters  7Z,  ui  and  x.  Thus  this  equation 
may  be  integrated  in  order  to  calculate  the  amplitude  : 


(3.1) 


With  the  same  meaning  as  used  in  the  boundary  layer,  see  [18],  the  factor  n  gives  the  logarithm 
of  the  eigenmode  amplitude.  In  equation  (3.1),  Xq  is  the  neutral  position.  If  it  exists,  it  is  the 
point  for  given  values  of  (7 Z,  to)  which  separates  upstream  a  stable  region  and  downstream  an 
unstable  region.  Generally,  Xq  depends  in  fact  on  (lZ,u>),  so  that  we  may  write  Xq  =  xo(1Z,u>).  In 
equation  (3.1)  Aq  is  an  amplitude,  which  corresponds  to  the  amplitude  at  the  abscissa  xo-  It  can 
be  noted  that  the  eigenmode  is  solution  of  a  homogeneous  problem,  so  that  its  definition  is  not 
univocal.  For  example,  with  the  equation  (2.9),  if  ip  satisfies  the  equation,  A %jj,  with  A  any  complex 
constant,  also  satisfies  the  equation.  For  that  reason,  the  amplitude  Aq  cannot  be  determined 
within  the  linear  stability  theory.  On  the  other  hand,  the  n  factor  is  intrinsic. 

3.2.2  Plane  case 

To  be  comparable  to  the  experimental  results  obtained  with  VECLA,  we  only  focus  on  the  amplified 
varicose  mode  in  this  paragraph.  To  give  an  example,  the  injection  Reynolds  number  TZ  is  fixed  : 
7 Z  =  4000.  Then,  the  growth  rate  —  oii  is  computed  for  different  frequencies,  and  for  different 
values  of  x.  Finally,  for  each  considered  frequency,  the  growth  rate  is  integrated  leading  to  the 
n  factor.  The  result  is  given  in  figure  3.4  in  a  diagram  (x,w).  Some  important  feature  may  be 
deduced  from  this  figure.  First  there  is  a  critical  abscissa,  which  is  close  to  5  for  the  considered 
Reynolds  number1.  It  means  that  upstream  the  abscissa  x  =  5,  the  Taylor  flow  is  stable,  possible 
eigenmodes  are  damped  in  this  region.  On  the  other  hand,  downstream  the  critical  abscissa,  i.e. 
for  x  >  5,  there  is  a  range  of  frequencies  for  which  the  instability  modes  are  amplified.  Moreover,  it 
can  be  noted  that  this  range  of  “dangerous”  frequencies  increases  with  x.  However,  low  frequencies 
(less  than  approximately  14)  are  never  amplified,  whereas  large  frequencies  become  amplified  but 
for  more  and  more  large  values  of  x. 

The  upstream  curve  corresponds  to  n  =  0,  i.e.  to  the  neutral  (or  marginal)  curve  :  the 
location  of  xq{TZ1u>)  with  the  notations  of  (3.1).  Without  any  additional  information,  it  is  usually 
assumed  that  the  amplitude  A0  is  a  constant  at  the  abscissa  x0(TZ,  u>).  In  fact  this  point  is  related  to 
the  receptivity  process  which  describes  the  physical  mechanism  by  which  eigenmodes  appear  from 
the  general  ambiant  turbulent  noise  or  from  any  other  types  of  excitations.  If  this  ambiant  noise 
does  not  contain  any  favoured  frequencies,  at  least  in  the  range  of  possible  amplified  frequencies, 
this  uniform  repartition  of  initial  amplitude  Aq  seems  natural  as  a  first  approximation.  In  this  case, 
the  n  factors  give  exactly  the  amplitude  of  the  eigenmodes  up  to  a  scale  factor.  Thus,  figure  3.4 

1In  fact  the  exact  value  of  the  critical  abscissa  must  be  carefully  considered.  This  is  due  to  the  non  parallel 
effects  which  strongly  increase  for  small  values  of  x,  see  section  4.1. 
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Figure  3.4:  Iso-n  factors  in  a  (x,u>)  diagram,  TZ  =  4000 


also  shows  that  the  frequency  which  corresponds  to  the  largest  eigenmode  at  a  given  distance  from 
the  front  wall  is  slightly  (and  quasi  linearly)  increasing  with  this  distance. 

Finally  this  figure  quantifies  the  amplification,  for  example,  at  x  =  20,  the  maximum  of  the 
n  factor  is  close  to  7,  it  means  that  the  perturbation  is  e7  ~  1000  times  larger  than  the  perturbation 
at  x  =  5. 

3.2.3  Axisymmetric  case. 

For  reasonable  values  of  the  injection  Reynolds  number,  the  frequency  and  the  values  of  x,  only 
one  mode  (for  each  value  of  the  azimuthal  wave  number)  has  been  found  to  become  amplified  from 
some  abscissa.  The  other  modes  remain  damped.  We  obviously  focus  on  this  possibly  amplified 
mode. 

Similarly  to  figure  3.4,  figure  3.5  gives  for  m  =  0  (axisymmetric  modes)  and  for  the  same 
Reynolds  number  as  before  1Z  =  4000,  the  values  of  the  computed  n  factor  in  a  diagram  (a:,  w).  The 
general  shape  is  similar  to  the  one  obtained  in  the  plane  case.  However  some  important  differences 
exist,  that  must  be  noted.  First  the  critical  abscissa  is  located  upstream  in  the  axisymmetric  case. 
Furthermore,  the  n  factors  seems  to  be  larger  in  this  geometry  than  in  the  plane  one.  For  example 
at  x  =  15,  the  n  factor  is  close  to  11  in  the  axisymmetric  geometry,  it  is  only  close  to  3  in  the 
plane  geometry.  Taken  into  account  that  the  growth  of  the  perturbation  follows  an  exponential 
behaviour,  this  means  that  the  axisymmetric  Taylor  flow  is  much  more  unstable  than  the  plane 
Taylor  flow.  This  difference  may  be  related  to  the  difference  observed  previously  in  figure  1.6 
between  the  streamlines  in  the  two  cases.  Indeed,  one  important  question  is  to  determine  the 
origin  of  the  present  intrinsic  instability.  Even  if  it  is  not  really  demonstrated,  it  may  be  suggested 
that  the  instability  comes  from  the  strong  curvature  of  the  streamline  close  to  the  porous  wall. 
This  could  be  at  least  coherent  with  the  following  observations. 
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Figure  3.5:  Iso-n  factor,  m  =  0,  1Z  =  4000 


1.  The  maximum  of  the  fluctuating  velocity  is  close  to  the  porous  wall. 

2.  The  Taylor  flow  is  stable  upstream  (where  the  curvature  of  the  streamlines  may  be  not  large 
enough) 

3.  If  the  non  parallel  terms  in  V  are  suppressed  in  the  stability  equations,  such  as  in  (2.8),  the 
basic  flow  is  found  to  be  always  stable,  see  [19]. 

4.  Finally,  on  one  hand,  the  curvature  is  much  more  stronger  in  the  axisymmetric  case  as  in  the 
plane  one,  as  observed  in  figure  1.6  and,  on  the  other  hand,  the  growth  rates  are  much  more 
larger  in  the  axisymmetric  geometry. 

The  last  difference  between  figures  3.4  and  3.5  concerns  the  frequency.  It  seems  that  the  amplified 
frequencies  are  greater  in  the  axisymmetric  case  than  in  the  plane  one.  For  example,  the  range  of 
amplified  frequencies  for  x  =  8  is  close  tou£  [30, 140]  whereas  it  is  only  u>  £  [15,40]  for  the  plane 
Taylor  flow. 

The  results  given  above  correspond  to  the  two-dimensional  (noted  2D)  mode,  that  is  for  the 
azimuthal  wave  number  m  =  0.  An  important  point  concerns  the  nature  of  the  most  amplified 
mode.  A  well-known  result  is  the  Squire’s  theorem,  which  expresses  that  for  a  2D  strictly  parallel 
basic  flow,  the  first  instability  occurs  for  a  2D  perturbation.  However,  it  does  not  mean  that  a 
two-dimensional  perturbation  is  always  more  amplified  than  any  three-dimensional  perturbation, 
see  [20]  for  more  information.  However  for  the  Taylor  flow  this  theorem  cannot  be  applied  even 
in  the  plane  case,  considering  that  the  stability  equations  contain  some  nonparallel  terms  which 
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makes  unapplicable  the  mentioned  theorem.  Some  direct  computations  in  the  plane  case  show 
however  that  the  most  amplified  modes  are  approximately  two-dimensional,  these  are  the  varicose 
and  the  sinuous  modes.  For  the  axisymmetric  geometry,  figure  3.6  shows  the  neutral  curve  and  the 


x 


Figure  3.6:  Marginal  curves  for  m  =  0, 1,  2,  3,  1Z  =  4000 


location  of  the  maximum  of  the  n  factor  in  the  ( x ,  to)  diagram  for  the  four  first  positive  azimuthal 
wave  numbers.  The  first  observation  is  that  the  results  for  to  =  0  and  to  =  1  are  quasi  identical. 
This  conclusion  exactly  corresponds  to  the  similarity  between  the  varicose  and  the  sinuous  modes 
previously  observed  for  the  plane  geometry.  On  the  other  hand,  larger  azimuthal  wave  numbers 
seem  to  be  less  and  less  amplified,  except  for  the  high  frequencies,  for  which  the  differences  seem 
to  be  small.  However  the  critical  abscissa  for  m  =  3  is  close  to  the  one  of  m  =  0,  for  example 
it  is  clearly  upstream  the  location  of  the  n  =  1  curve  in  figure  3.5.  This  means  that  if  all  the 
eigenmodes  (characterised  by  (m,aj))  start  with  the  same  initial  amplitude  Aq  at  their  neutral 
curve,  the  amplitude  of  the  m  =  0  mode  at  the  critical  abscissa  of  m  =  3  is  only  a  little  larger 
than  the  amplitude  of  the  m  =  1,2,3  modes.  Therefore,  it  is  clear  that  at  a  given  x  position, 
x  =  7  for  example,  the  perturbation  should  be  the  superposition  of  several  modes  corresponding 
to  different  values  of  to,  so  that  for  any  hot  wire  measurement  it  should  be  very  difficult  to  isolate 
the  contributions  of  the  different  modes. 

To  conclude  this  overview  of  the  typical  stability  results,  it  remains  to  analyse  the  real  part 
of  the  complex  wave  number  a.  It  is  then  more  suitable  to  work  with  the  wavelength  A  which 
is  defined  by  :  A  =  2-7r/ar.  Figure  3.7  gives  the  iso-A  values,  always  in  a  (x,u>)  diagram  for  the 
axisymmetric  geometry  with  TZ  =  4000  and  m  =  0.  It  is  interesting  to  note  that  at  least  for 
values  of  x  less  than  15,  the  wavelength  of  the  mode  which  corresponds  to  the  largest  amplitude 
is  close  to  1.  With  a  dimensional  point  of  view,  this  means  that  the  largest  mode  presents  a 
wavelength  roughly  equal  to  the  radius  of  the  duct.  In  particular  the  streamwise  evolution  of  the 
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Figure  3.7:  Iso- wavelength  for  m  =  0,  1Z  =  4000 


eigenmodes  occurs  on  a  scale  which  is  much  shorter  than  the  total  length  of  the  duct  which  scales 
the  undesirable  cavity  modes.  Due  to  this  great  difference  of  streamwise  extend,  it  may  be  guessed 
that  the  instability  alone  is  not  very  efficient  as  a  possible  exciting  source  of  the  cavity  modes. 


3.3  Influence  of  the  Reynolds  number 

Concerning  the  dispersion  relation,  the  influence  of  the  different  parameters  x,  to  and  the  azimuthal 
wave  number  for  the  axisymmetric  geometry  have  been  studied.  It  remains  the  influence  of  the 
injection  Reynolds  number.  For  example,  the  plane  geometry  is  considered,  with  fixed  x  and  ui 
values  :  x  =  10,  to  =  31.6  (which  have  been  chosen  in  the  spectrum  analysis  of  the  linearised 
operator,  see  figure  3.1).  The  result  is  plotted  in  figure  3.8  for  the  streamwise  wave  number  ar 
represented  in  dashed  line  and  for  the  growth  rate  —  a*  represented  in  full  line.  In  this  figure,  the 
basic  flow  is  always  the  one  given  by  the  analytical  Taylor  form  (1.8)  and  only  for  the  stability  the 
Reynolds  has  been  varied  from  50  up  to  105.  It  is  clear  that  there  is  an  asymptotic  behaviour  when 
the  Reynolds  number  increases.  This  means  that  the  basic  mechanism  of  this  instability  is  mainly 
inviscid.  Up  to  now,  a  clear  physical  explanation  “with  the  hands”  of  the  origin  of  the  instability 
has  not  be  found.  As  explained  in  the  previous  section,  the  instability  appears  to  be  related  to 
the  strong  curvature  of  the  streamlines  of  the  mean  flow  close  to  the  injecting  wall,  which  is  little 
affected  by  the  viscosity. 

However  the  wave  number  seems  to  tend  to  the  inviscid  value  faster  than  the  growth  rate. 
Thus,  especially  for  the  estimation  of  the  n  factor  in  practical  case  for  which  the  Reynolds  number 
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Figure  3.8:  Influence  of  the  Reynolds  number  on  the  complex  streamwise  wave  number  a,  plane 
case,  x  =  10,  u>  =  31.6,  varicose  mode.  The  growth  rate  —  oti  is  in  full  line,  the  tenth  of  the  wave 
number  oy/10  in  dashed  line 


is  of  order  1000,  it  is  necessary  to  take  into  account  the  viscous  effects,  the  difference  could  be  of 
about  10%  for  inviscid  computations.  On  the  other  hand,  it  is  important  to  emphasize  that  for 
the  mean  flow,  it  is  really  not  necessary  to  use  the  viscous  mean  flow  given  by  (1.7)  for  Reynolds 
numbers  greater  than  100,  the  Taylor  solution  is  accurate  enough. 

The  same  feature  can  be  observed  in  the  axisymmetric  geometry,  leading  to  the  same  con¬ 
clusions. 


3.4  Comparisons  with  the  experiment 

3.4.1  Preparation  of  the  results  for  the  comparison 

Hot  wire  measurements  give  access  to  the  instantaneous  velocity  at  the  position  where  the  probe 
has  been  placed.  Then,  using  different  possible  treatments  of  the  signal,  the  latter  is  decomposed 
in  two  parts  :  the  mean  value  and  the  fluctuating  one  upon  which  a  Fourier  transform  is  applied  in 
order  to  get  the  spectral  dependence  of  the  fluctuating  velocity  or  more  often  the  power  spectral 
density  (DSP). 

This  spectral  representation  may  be  compared  to  the  theoretical  predictions.  For  each  given 
frequency,  the  n  factor  is  computed  by  integrating  in  x  the  amplification  growth  rate,  so  that  the 
amplitude  of  the  fluctuation  is  simply  theoretically  given  by  HoGxp(n(x,w)),  with  A$  the  initial 
amplitude  coming  from  the  receptivity  conditions.  In  the  following  results,  this  constant  A$  is 
assumed  to  be  independent  of  ui,  but  to  be  a  function  of  the  experimental  conditions  (values  of  the 
injection  velocity,  of  the  height  of  the  duct,  of  the  porosity  of  the  injecting  walls,  ...).  Obviously,  Aq 
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is  independent  of  x.  Consequently  for  given  experimental  conditions,  the  constant  Aq  is  adjusted 
once,  in  order  to  fit  for  one  frequency  and  one  value  of  x  the  measured  amplitude.  In  case  of 
comparisons  including  y  variations,  the  same  constant  Aq  must  be  considered  of  course. 

3.4.2  Effects  of  the  dimensional  height  of  the  duct  and  injection  velocity 

It  is  interesting  to  test  first  the  scaling  effects  by  varying  experimentally  the  injection  velocity  and 
the  height.  From  the  theoretical  point  of  view,  as  the  computation  of  the  n  factor  in  the  (x,oj) 
plane  is  performed  with  dimensionless  quantities,  changing  the  injection  velocity  and  the  height 
only  modifies  the  Reynolds  number.  Comparisons  with  measurements  obtained  with  the  VECLA 
facility  are  reported  in  figure  3.9.  Four  comparisons  are  given,  they  correspond  to  four  values  of  the 
pair  (injection  velocity  Vjnj,  height  h).  The  measured  values  are  plotted  in  full  line,  the  theoretical 
predictions  by  the  open  circles.  In  the  four  cases,  the  frequency  is  plotted  between  0  and  2000  Hz. 
However  the  scale  of  the  vertical  axis,  giving  the  amplitude  of  the  fluctuation,  differs  from  each 
other.  The  adjusted  initial  amplitude  Aq  is  also  given  for  each  comparison  in  the  legend,  it  seems 
that  the  constant  Aq  increases  with  the  injection  velocity  and  decreases  with  the  height  of  the 
duct. 

The  overall  comparison  is  quite  satisfying  (good  prediction  of  the  amplified  frequencies  and 
quite  good  shape  of  the  spectral  dependence) ,  even  if  there  are  some  differences  which  are  mainly 
visible  for  h  =  10  mm.  It  may  be  remarked  however  that  for  this  height  (the  smallest  one  for  the 
VECLA  facility),  the  mean  flow  does  not  exactly  coincide  with  the  theoretical  Taylor  flow  (or  the 
viscous  solution),  it  seems  that  small  three-dimensional  structures  slightly  modify  the  streamwise 
component  of  the  mean  velocity  in  comparison  with  the  expected  one.  However  the  broadness 
of  the  amplified  frequencies  peak  as  well  as  the  frequency  associated  to  this  peak  are  in  good 
accordance  with  the  linear  stability  results. 

A  similar  comparison  has  been  done  for  the  axisymmetric  geometry  with  the  available  ex¬ 
perimental  results,  see  [5].  However  in  this  case,  the  experiment  have  been  carried  out  in  1990, 
that  is  before  the  knowledge  of  the  intrinsic  instability  described  in  the  present  document.  Among 
other  comparisons,  which  can  be  found  in  [16],  figure  3.10  shows  the  power  spectral  density  of 
the  fluctuating  axial  velocity  for  three  values  of  the  injection  velocity.  The  latter  is  represented 
by  Mw,  the  wall  injection  Mach  number.  The  left  figures  (experiment)  are  directly  scanned  from 
the  publication,  whereas  the  right  ones  give  the  theoretical  results  obtained  for  the  four  first  posi¬ 
tive  azimuthal  wave  numbers.  Three  different  initial  amplitudes  Aq  have  been  chosen  for  the  three 
different  injection  velocities.  As  before,  the  range  of  the  amplified  frequencies  seem  to  agree  quite 
well  with  the  experimental  results,  even  if  there  are  some  discrepancies  which  are  significantly 
larger  than  the  ones  observed  with  VECLA  in  the  plane  geometry. 

3.4.3  Streamwise  amplification 

The  previous  comparisons  proved  the  good  agreement  between  the  experimental  results  and  the 
theoretical  ones  in  terms  of  range  of  amplified  frequencies.  The  goal  of  this  section  is  to  analyse 
carefully  the  instability  process  itself,  that  is  the  streamwise  amplification,  which  is  theoretically 
of  exponential  type.  The  fluctuating  velocity  at  different  values  of  x  must  be  then  compared  in 
the  same  experimental  condition  (in  order  to  keep  the  same  initial  amplitude  Aq  for  the  different 
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Without  nozzle,  porosity  =  18  pun,  y0=  1  mm 


Vi„j=1.36m/s,  h=10mm 


Vtoj=1.36m/s,  h=20mm 

*10’ 


Vinj=  1 . 70m/s,  h=l  Omm 

*10’ 


V^=1.70m/s,  h=20mm 

*10’ 


Figure  3.9:  Comparisons  between  LST  results  and  experimental  ones  with  the  VECLA  facility 
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Mw  =  0.0036 


Mrt  =  0.0027 


Mn.  =  0.0018 


Figure  3.10:  Comparisons  between  LST  results  and  experimental  ones  obtained  by  Brown  and 
co-workers 
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values  of  x). 

A  first  comparison  is  given  in  figure  3.11.  The  height  of  the  duct  is  30  mm.  The  experimental 
results  (lines)  are  given  in  six  increasing  abscissa,  whereas  the  theoretical  values  (symbols)  are 
reported  for  only  the  four  largest  values  of  x.  The  two  first  values  are  indeed  upstream  the  critical 
abscissa.  The  comparison  gives  an  excellent  agreement  between  the  experiment  and  the  theory 
demonstrating  undoubtedly  on  one  hand  the  real  existence  of  intrinsic  instability  for  this  type  of 
flow  and  on  the  other  hand  the  reliability  of  the  present  linear  stability  theory.  To  emphasize  this 
comparison,  let  us  mention  :  the  range  of  amplified  frequencies  (as  shown  above),  the  streamwise 
amplification  (as  previously  said,  the  same  constant  Aq  =  1/6000  has  been  used  for  the  different 
abscissa).  Furthermore,  the  increasing  shift  of  the  amplified  frequencies  with  respect  to  x  is  also 
clearly  visible  in  this  figure.  Finally,  the  spectrum  measured  at  x  =  570  mm,  the  most  downstream 


*10' 2  VECLA,  height  =  30  mm,  Vinj  =  1.36  m/s 


Figure  3.11:  Comparisons  between  LST  results  and  experimental  ones  with  the  VECLA  facility 


position,  exhibits  an  interesting  (weakly)  nonlinear  phenomenon,  that  is  the  increase  of  2w  and  0 u> 
(steady)  modes.  These  first  nonlinear  mechanisms  are  quite  well  understood  in  the  framework  of 
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the  boundary  layer  stability,  see  [21].  The  basic  reason  comes  from  the  linear  form 

uv  =  -(«exp(iw)  +  u*  exp 

of  the  physical  fluctuation  (where  z*  represents  the  complex  conjugate  of  z).  The  quadratic  terms  in 
the  Navier-Stokes  equations  are  formally  expressed  by  the  product  u^u^.  The  product  is  converted 
into  a  sum  of  the  exponential  terms,  leading  thus  to  the  frequencies  2 u>  and  0.  And,  as  the  linear 
step  does  not  amplify  only  one  mode,  but  a  rather  broad  range,  the  quadratic  interactions  lead  also 
to  amplify  a  certain  range  around  zero  and  another  one  corresponding  to  twice  the  linear  range. 
Similar  comparisons  can  be  also  tested  with  the  axiymmetric  duct.  In  the  present  document, 


Figure  3.12:  Comparisons  between  LST  results  and  experimental  ones  obtained  with  the  facility 
of  the  VKI.  The  full  lines  correspond  to  the  experimental  results,  each  of  them  at  one  value  of  x. 
They  are  ordered  regularly  :  the  large  the  FFT  amplitude,  the  larger  the  value  of  x. 


two  different  comparisons  are  shown,  the  experimental  results  coming  from  two  different  facilities. 
The  first  one  is  located  in  VKI  and  the  results  have  been  obtained  by  J.  Anthoine,  see  [6].  As 
before,  the  axial  fluctuating  velocity  spectrum  is  compared  between  both  approaches,  the  results 
are  given  in  figure  3.12.  The  injection  velocity  is  0.88  m/s,  the  radius  of  the  duct  38  mm.  The 
spectra  are  given  for  three  values  of  x,  between  7.4  and  10  in  terms  of  dimensionless  distances  from 
the  front  wall.  In  comparison  with  the  result  given  for  the  plane  geometry,  it  appears  immediately 
that  the  measured  velocity  is  much  less  smooth  than  in  figure  3.11. 
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Similar  behaviour  have  been  recently  obtained  by  using  the  facility  VALDO,  see  [22],  The 
results  are  reported  in  figure  3.13.  The  radius  is  30  mm,  the  injection  velocity  1.05  m/s,  the  chosen 
abscissa  are  thus  close  to  10.5  in  dimensionless  value. 

From  both  figures,  it  may  be  concluded  that  there  is  a  range  of  amplified  frequencies  which 
correspond  roughly  to  the  theoretical  predictions.  It  seems  also  that  there  is  a  streamwise  amplifica¬ 
tion  but  some  discrepancies  are  clearly  present  between  the  theoretical  results  and  the  experimental 
ones.  It  is  not  easy  to  explain  the  reason  of  this  noisy  signal  obtained  in  the  axisymmetric  facil- 
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Figure  3.13:  Comparisons  between  LST  results  and  experimental  ones  obtained  with  the  VALDO 
facility.  Each  curve  corresponds  to  one  of  the  four  measurements,  the  same  is  true  for  the  symbols 
giving  the  theoretical  predictions 


ity.  This  may  be  due  to  the  existence  of  several  modes  (associated  to  different  azimuthal  wave 
numbers),  but  this  may  be  also  due  to  the  measurement  anemometer.  In  the  plane  case,  the  hot 
wire  is  parallel  to  the  porous  wall  whereas  in  the  other  geometry,  some  corners  effects  may  occur. 
As  the  probe  is  placed  very  close  to  the  porous  boundary,  there  are  possible  interactions  between 
the  mean  flow  and  the  corners  of  the  probe  which  carry  the  hot  wire  itself.  Anyway,  it  may  be 
noted  that  the  streamwise  amplification  is  very  large.  For  example,  the  result  shown  in  figure  3.13 
exhibits  a  huge  amplification  (factor  about  20  in  the  experiment  over  a  distance  of  only  1.5  cm). 
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Chapter  4 

In(tro)spection  of  the  used 
assumptions 


Different  assumptions  have  been  made,  their  influences  must  be  inspected  in  order  to  determine 
the  generality  of  the  results  presented  in  the  previous  chapter. 


4.1  Non  parallel  effects 

4.1.1  Definitions 

As  explained  before,  the  basic  flow  explicitly  depends  on  the  streamwise  coordinate  x  so  that 
the  exponential  form  exp(ia:r)  used  for  the  perturbation  is  not  theoretically  justified.  For  a  two- 
dimensional  basic  flow,  a  dependence  of  the  streamwise  velocity  coordinate  with  respect  to  x  is 
related  to  a  non  zero  transverse  velocity  V  through  the  continuity  equation.  This  is  the  reason 
why  this  x  dependence  is  called  a  non  parallel  effect.  In  the  following,  the  basic  flow  is  assumed 
to  be  dependent  on  x  (even  weakly). 

For  clarity,  let  us  call  the  parallel  approach  the  one  when  only  the  U  is  kept  in  the  stability 
equations  (leading  to  exactly  the  Orr-Sommerfeld  equation)  and  the  nonparallel  approach  the  one 
for  which  all  terms  related  to  the  basic  flow  are  kept  in  the  stability  equations.  Then,  if  the 
perturbation  has  the  form  of  a  normal  mode  for  a  parallel  basic  flow  (or  for  a  flow  of  which  the 
non  parallel  terms  are  neglected),  the  corresponding  approach  will  be  noted  by  OSE  (for  Orr- 
Sommerfeld  Equation).  On  the  other  hand,  the  (theoretically  not  justified)  use  of  the  normal 
mode  for  a  non  parallel  basic  flow  will  be  noted  by  NNP  (for  Normal  Non  Parallel  approach).  All 
the  results  shown  in  the  previous  chapter  have  been  obtained  by  using  the  NNP  approach. 

Now  there  are  some  questions  :  is  the  Taylor  flow  weakly  non  parallel  (so  that  OSE  is  a 
more  or  less  accurate  approximation)  ?  Is  it  possible  to  justify  the  use  of  NNP  ?  Is  it  possible 
to  perform  a  consistent  and  accurate  non  parallel  stability  analysis  ?  The  following  part  of  the 
present  chapter  is  aimed  to  provide  some  answers  to  that  questions. 
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4.1.2  Use  of  a  parallel  approximation  ? 


As  explained  above,  the  non  parallel  effects  are  connected  with  the  magnitude  of  V  with  respect 
to  U.  Using  the  Taylor  analytical  form  (1.8),  the  ratio  V /U  may  be  readily  calculated  : 


V  2  I  sin  ^  I  2  n\y\ 

=  —  J - !-!  =  —  tan 

U  7 TX  COS  ~Y  TTX  2 


This  suggests  that  this  ratio  varies  from  0  (at  the  symmetry  axis)  to  oo  at  the  porous  walls. 
However  if  the  dependence  on  y  is  suppressed  by  a  norm  such  as  the  maximum  with  respect  to  y, 
the  ratio  is  simplified  into  V/U  =  2/itx,  which  indicates  that,  according  to  this  choice  of  norm, 
the  non  parallel  effects  decrease  with  respect  to  x  (the  ratio  is  smaller  than  0.05  for  x  >  15  for 
example).  From  these  simple  two  remarks,  it  can  be  concluded  that  the  Taylor  flow  may  appear 
more  and  more  parallel  for  large  values  of  x  but,  reciprocally,  at  any  values  of  x  it  exists  a  small 
region  close  to  the  injecting  walls  where  the  basic  flow  is  strongly  non  parallel. 

In  order  to  illustrate  the  non  parallel  effects,  a  first  possibility  consists  in  solving  the  OSE 
with  the  streamwise  velocity  component  U  being  the  one  given  by  Taylor,  even  if  the  mean  flow 
given  by  the  velocity  (7ra;/2  cos(7n//2),  0)  obviously  does  not  satisfy  the  steady  equations.  Figure  4.1 
gives  for  the  axisymmetric  Taylor  flow,  with  7Z  =  4500,  q  =  0  and  ui  =  80,  the  evolution  with  respect 


Figure  4.1:  Comparison  of  OSE,  NNP  and  PSE  results  for  the  axisymmetric  Taylor  flow,  1Z  =  4500, 
q  =  0  and  u  =  80.  The  wavenumber  ar  is  plotted  in  the  left  hand  side,  the  growth  rate  at  in  the 
right  hand  side 


to  x  of  the  complex  number  a  solved  by  the  OSE,  the  NNP  and  the  PSE  approaches  (see  below 
for  some  explanations  on  the  PSE  method).  It  appears  that  on  one  hand  NNP  and  PSE  results 
are  very  similar,  whereas  OSE  strongly  differs  from  the  two  others,  for  example  the  growth  rate 
remains  always  positive.  Similar  behaviour  has  been  found  in  plane  geometry,  see  [19]. 

This  strong  non  parallel  effect  is  also  confirmed  by  the  whole  spectrum  obtained  with  the 
OSE  linearised  operator  with  the  same  values  of  the  parameters  are  for  figure  3.1.  The  results 
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are  given  in  figure  4.2.  Once  more,  the  real  part  is  completely  different  and  no  amplified  mode 
5 1 - i - 1 - 1 - 1 - 


Figure  4.2:  Spectrum  of  the  artificially  parallel  Taylor  flow,  1Z  =  1000,  x  =  10,  a  =  4 


has  been  found.  Let  us  recall  that  the  experiment  clearly  shows  amplified  modes  which  are  in 
relative  good  agreement  with  the  NNP  results.  The  parallel  OSE  approach  does  not  seem  to  be 
an  approximation  of  the  stability  results  for  the  present  basic  flow. 

4.1.3  What  is  the  matter  with  the  normal  non  parallel  approach  ? 

As  firstly  remarked  by  comparing  results  between  those  of  Varapaev  and  Yagodkin  [7] ,  those  of  Lee 
and  Beddini  [23]  and  ours  [19],  there  is  a  problem  somewhere.  Finally  the  explanation  has  been 
found  by  Griffond  [24],  the  NNP  approach  is  not  consistent,  because  the  stability  problem  depends 
on  the  formulation  which  is  used.  In  the  present  course,  two  formulations  have  been  used  in  the 
plane  case  :  with  the  primitive  variables  ( u,v,p )  see  system  (2.8)  and  with  the  stream  function  ip, 
see  equation  (2.9).  The  matter  is  to  determine  if  both  formulations  are  strictly  equivalent.  Before 
introducing  the  normal  mode,  both  formulations  (1.1)  and  (1.3)  are  obviously  equivalent  as  well  as 
the  linearised  equations  deduced  from  each  of  them.  But,  as  the  normal  mode  is  not  and  cannot  be 
justified,  the  question  of  the  equivalence  between  system  (2.8)  and  equation  (2.9)  remains  posed. 

The  answer  given  in  [24]  is  the  following  one.  From  system  (2.8),  the  amplitude  function  of 
the  pressure  p  may  be  eliminated  by  multiplying  the  third  equation  of  that  system  by  ia,  differenti¬ 
ating  the  second  one  (with  respect  to  y)  and  finally  by  subtracting  both  obtained  equations.  Then 
in  this  new  equation,  the  amplitude  functions  u  and  v  may  be  replaced  by  respectively  —iaip  and 
Dip  where  ip  is  an  amplitude  function  for  the  stream  function  and  D  the  differentiation  operator 
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with  respect  to  y.  If  we  proceed  in  this  way,  the  following  equation  is  obtained 

/  1  \  „  >)2r y  pp. fj  r  arr  A 

(  -iu)  +  iaU  +  VD  -  —  {D2  -  a2)  )  ( D 2  -  a2)ip  — — — .Eh/S  -  +  —a2^  =  0  (4.1) 

\  7c  )  oyz  oyz  ox 

which  is  NOT  equation  (2.9)  as  it  ought  to  be  in  a  consistent  theory,  the  last  term  of  (4.1)  is  not 
present  in  (2.9).  Let  us  just  recall  that  for  a  strictly  parallel  basic  flow,  both  formulations  are 
strictly  equivalent,  they  lead  both  to  the  Orr-Sommerfeld  equation. 

The  NNP  approach  is  consequently  not  consistent,  its  results  may  depend  on  the  formulation. 
To  illustrate  that  point,  comparisons  between  the  two  formulations  have  been  performed,  they  are 
illustrated  in  figure  4.3.  In  addition,  the  published  results  of  [7]  and  [23]  have  been  indicated  in 


Figure  4.3:  Comparisons  between  diverse  published  results  and  our  results  illustrating  the  incon¬ 
sistency  of  the  normal  non  parallel  approach.  Using  the  stream  function  formulation,  results  for 
the  neutral  curve  for  the  sinuous  modes  in  [7,  23]  (circles)  and  for  the  varicose  modes  in  [23]  coin¬ 
cide  with  ours  (full  lines).  They  are  significantly  different  from  ours  obtained  with  the  primitive 
variable  formulation  for  the  varicose  mode  (dashed  line)  and  for  the  sinuous  mode  (dotted  line). 
Self-similar  solution  is  used  for  the  plane  basic  flow  with  1Z  =  100. 


symbols.  The  neutral  curve  has  been  calculated  for  1Z  =  100  for  both  sinuous  and  varicose  modes 
and  with  the  two  formulations.  It  is  represented  in  the  ( Rc  ,  ar)  diagram,  where  Rc  represents  the 
Reynolds  number  based  on  the  axial  velocity  (which  is  linearly  dependent  on  x) .  All  the  published 
results  seem  to  be  accurately  calculated,  the  differences  between  the  full  lines  and  the  dashed  or 
the  dotted  ones  are  intrinsic,  they  are  due  to  the  choice  of  the  normal  mode  which  is  inconsistent 
with  the  present  basic  flow. 
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4.1.4  Is  there  any  consistent  and  accurate  non  parallel  approach  ? 


From  a  practical  point  of  view,  the  inconsistent  approach  seems  however  to  work  with  enough 
accuracy.  For  large  Reynolds  number  as  well  as  for  large  values  of  x,  the  differences  between  the 
two  approaches  discussed  in  the  previous  paragraph  are  quasi  negligible.  Furthermore,  as  shown 
in  the  previous  chapter,  the  agreement  with  the  experimental  results  is  very  good.  Even  for  the 
axisymmetric  geometry,  for  which  the  growth  rates  are  larger  and  the  involving  x  values  smaller,  it 
may  be  guessed  that  for  practical  prediction  the  NNP  approach  provides  a  satisfactory  estimation 
of  the  instability. 

On  the  other  hand,  the  present  solution  is  clearly  wrong  for  the  small  values  oi  x.  In 
particular  the  location  of  the  neutral  curve  must  be  considered  with  care.  Thus  an  alternative 
approach  for  the  mathematical  form  of  the  perturbation,  which  would  be  more  consistent  and 
more  accurate  with  respect  to  the  non  parallel  characteristics  of  the  basic  flow,  would  be  useful. 

Some  approaches  already  exist  which  are  first  consistent  and  whose  goal  secondly  is  to  deal 
with  the  non  parallel  effects.  The  clearest  one  is  based  on  the  multiple  scale  analysis,  see  [25] 
or  [26]  for  a  general  presentation  and  see  [27]  for  an  example  of  use  in  a  non  parallel  stability 
analysis  (in  the  framework  of  the  boundary  layer  stability).  A  more  recent  approach,  called  PSE 
(for  Parabolized  Stability  Equations),  developed  by  Herbert  and  Bertolotti  [28],  also  deals  with 
the  non  parallel  effects.  This  approach  has  been  used  for  the  present  mean  flow,  the  results  are 
shown  in  figure  4.1.  They  prove  that  after  a  numerical  transient  in  x,  the  PSE  results  rapidly  agree 
with  the  NNP  ones.  Concerning  the  consistency,  the  multiple  scale  analysis  is  fully  consistent,  the 
PSE  approach  is  quasi  consistent,  see  [29]  for  more  details.  The  (theoretical)  problem  with  these 
approaches  is  that  both  consider  the  Orr-Sommerfeld  equation  as  a  first  order  solution,  whereas  it 
is  not  the  case  for  the  Taylor  basic  flow. 

Another  mathematical  form  for  the  perturbation  has  been  proposed  by  J.  Griffond,  see  [24]. 
The  starting  equation  is  the  linearised  form  of  the  Navier-Stokes  equations  expressed  with  the 
stream  function  (1.3).  Using  the  perturbation  technique  ip  =  ip  +  ip,  the  linearised  problem  is  : 


d_ 

at 


Aip  + 


dip  d 
dy  dx 


Aip  + 


dip 

dy 


dip  d  d  -dip 

—  — A  ip - Aip  — 

dx  dy  V  dy  V  dx 


— AA  ip 

n  s 


(4.2) 


with  the  following  boundary  conditions  : 


V*/  e  [-1,4-1]  ^(o,j/)  =  o 

dip ,  .  dip .  .  dil’ ,  .  dip ,  . 

Vz  >  0  —  (x,  -1)  =  0,  -^-(x,l)=Q  —  (x,  1)  =  0 


(4.3) 


As  shown  before,  a  possible  choice  for  one  stream  function  associated  to  the  Taylor  basic  flow 
is  ip  =  xsm.^Try/2).  Equation  (4.2)  represents  then  a  linear  equation  with  coefficients  which  are 
dependent  on  x  and  y  and  independent  of  t.  Thus  a  general  form  for  the  unknown  function  ip(x ,  y,  t) 
may  be  chosen  with  an  exponential  dependence  with  respect  to  t  (normal  mode  in  t)  : 


ip{x,y,t)  =  e  lu}tip(x,y) 


Equation  (4.2)  becomes  then  : 
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whereas  the  associated  boundary  equations  remain  the  same  as  (4.3),  with  ip  instead  of  ip.  A 
possible  solution  could  be  to  solve  directly  the  stability  problem  given  by  the  previous  homogeneous 
partial  differential  equation.  Such  problems  are  currently  investigated  by  some  researchers  (in  case 
of  the  stability  of  a  recirculating  zone,  a  vortex,  the  boundary  layer  in  the  vicinity  of  the  attachment 
line  etc.)  but  there  is  a  strong  difficulty  with  the  treatment  of  the  boundary  conditions.  The 
proposed  solution  is  to  search  ip  in  the  form  : 

tp  =  xxip(y) 


leading  hence  to  ordinary  differential  equations  (like  with  the  normal  mode)  and  avoiding  therefore 
the  difficulty  with  the  boundary  conditions.  The  idea  with  this  mathematical  form  is  to  define  a 
phase  velocity,  which  is  then  uniformly  accelerated  with  x,  as  is  the  basic  Taylor  flow,  whereas 
with  the  normal  mode,  the  phase  velocity  does  not  depend  on  x.  With  this  form,  equation  (4.2) 
writes  : 


0  =  -  —  DiPp-iuD2iP  + 


Xtt  ttv  Try  „ 

—  cos  — - sin  —.D 

2  2  2 


D2ip  - 


-ip 
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1  1 

"ad  K 


- D  —  iu>  - 

n 


7 r  ny  . ,  .  .  Try  _ 

—  cos  —  (A  —  2)  —  sin  —  .D 
2  2  2 


A(A  -  l)i/> 


(4.4) 


A(A  —  1)(A  —  2)(A  —  3)ip 


with  D  the  differentiation  operator  with  respect  to  y.  As  it  is  the  case  with  the  normal  mode,  this 
form  does  not  satisfy  exactly  the  equations,  so  that  A  is  not  really  independent  of  x.  By  using  the 
so-called  quasi-parallel  approach,  the  value  of  A  is  recalculated  at  each  considered  value  of  x. 

The  advantage  of  this  form  is  clearly  to  define  a  zero-order  approximation  which  is  asymp¬ 
totically  valid  for  large  values  of  x  : 


- DAip  —  iojD2ip  + 
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On  the  other  hand,  let  us  return  to  the  OSE  and  NNP  approaches 
exp^cra))  for  which,  as  proposed  by  J.  Griffond,  we  define  : 


7 r2  -A 

—iP)  =  0  (4.5) 

(obtained  with  the  term 


A  =  ixa 


Then,  the  OSE  approach  gives  : 
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whereas  the  NNP  approach  leads  to  : 
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ip  =  0 


These  equations  show  that  for  large  values  of  x,  the  leading  order  (4.5)  of  (4.4)  coincides  exactly 
with  the  leading  order  of  the  NNP  approach  and  NOT  with  the  OSE  approach.  The  inconsistent 
approach  remains  inconsistent  obviously  but  is  now  justified  (its  leading  order  terms)  for  large 
values  of  x.  Finally  some  comparisons  between  direct  simulations  (by  using  the  code  SIERRA 
developed  by  F.  Vuillot,  see  [30])  and  stability  computations  show  a  better  agreement  with  (4.4) 
than  with  the  NNP,  see  [24]. 
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4.2  Physical  assumptions 

In  addition  to  the  assumption  of  a  normal  mode,  there  are  several  other  aspects  which  have  been 
neglected  in  the  presentation  given  before.  However,  the  effects  of  : 

•  compressibility, 

•  real  geometry, 

•  presence  of  an  acoustic  mode, 

•  regression  of  the  boundary  wall 

have  been  found  to  be  nearly  completely  negligible.  The  intrinsic  instability  described  before  seems 
to  be  very  robust. 

The  last  point  related  to  the  linear  stability  theory,  which  appears  to  be  important  in 
the  computations,  see  [31],  concerns  the  effects  of  the  presence  of  particles,  especially  of  reactive 
particles.  A  first  attempt  in  that  direction  has  been  achieved  by  taking  into  account  non  reactive 
particles.  A  first  interesting  result  of  [32]  is  that  the  particles  may  destabilise  the  flow,  as  it 
is  illustrated  in  figure  4.4.  This  effect  may  be  attributed  to  the  appearing  augmentation  of  the 


Figure  4.4:  Evolution  of  the  growth  rate  a*  as  function  of  the  density  of  the  injected  particles  at 
the  wall.  The  injection  velocity  of  the  particles  is  the  same  as  the  one  of  the  fluid,  the  Stokes 
number  is  10-3,  and  x  =  10,  u  =  30. 


Reynolds  number  due  to  the  augmentation  of  the  density  thanks  to  the  particles. 
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In  addition  to  this  two-phase  flow  effect,  there  is  of  course  the  linear  assumption.  But 
nonlinear  mechanisms  induced  by  instability  are  neither  very  easy  nor  very  common.  This  could 
be  another  lecture  !  See  the  thesis  of  J.  Griffond  [9]  for  detailed  nonlinear  analysis. 
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Conclusion 


Undoubtedly,  an  intrinsic  instability  exists  in  the  flow  induced  by  wall  injection.  A  linear  stability 
theory  can  be  carried  out,  even  if  it  presents  some  difficulties  (mainly  from  the  theoretical  point 
of  view)  in  relation  with  the  non  parallel  characteristics  of  the  considered  mean  flow.  At  least  in 
the  plane  geometry,  the  experimental  results  confirm  the  predictions  based  on  the  linear  stability 
theory.  However  in  practice,  the  intrinsic  perturbation  cannot  grow  infinitely  in  x ,  as  it  could  be 
according  to  the  linear  stability  results.  From  flow  visualisations  performed  in  VECLA,  we  know 
that  large  structures  are  emitted  somewhere,  they  are  associated  to  fluctuations  of  large  amplitude 
and  are  thus  not  compatible  with  a  small  perturbation  approach  leading  to  the  linear  theory.  The 
perturbations  grow  in  x,  according  to  the  linear  instability  mechanism  and  then  become  fully  non 
linear.  Moreover,  for  some  configurations  the  experimental  results  of  VECLA  do  not  seem  to 
exhibit  any  instability  feature.  Among  these  results,  there  is  the  resonance  on  an  acoustic  mode 
of  large  amplitude  which  is  of  course  the  most  important  configuration  in  practice. 

Where  does  this  resonance  phenomenon  come  from  ?  Does  the  stability  play  any  role  ? 

In  order  to  try  to  give  some  answers  to  these  questions,  let  us  summarise  the  main  ideas 
(coming  from  the  experimental  results,  direct  simulations  and  nonlinear  stability  theories)  obtained 
at  the  end  of  the  PhD  of  J.  Griffond.  The  key  point  is  the  dimensional  value  of  the  abscissa  of  the 
motor  exit  section  x*.  It  seems  that  the  abscissa  related  to  the  instability  mechanisms  (location 
of  the  iso-n  factor  curves,  location  of  the  large  structure  emission)  are  mainly  constant  in  term 
of  dimensionless  values,  with  the  reference  scale  for  the  distance  being  the  radius  of  the  fluid  h. 
But  for  a  real  motor,  only  x*  remains  constant,  the  value  of  h  grows  continuously  thanks  to  the 
regression  of  the  injecting  wall  due  to  the  combustion. 

The  proposed  scenario  is  then  the  following  one  : 

•  At  the  beginning,  h  is  small.  In  that  case,  the  flow  becomes  turbulent  somewhere  in  the  duct, 
the  large  structures  which  are  emitted  upstream  in  the  laminar  zone  are  dissipated  by  the 
turbulence,  so  that  there  is  no  coherent  (i.e.  quasi  periodic  in  time)  structure  which  passes 
through  the  exit  section.  Nothing  happens  with  respect  to  the  cavity  modes. 

•  Later,  h  becomes  larger,  the  transition  to  turbulence  moves  downstream  and  at  a  given  time,  the 
turbulence  zone  is  not  sufficient  to  dissipate  enough  the  large  emitted  structures.  The  frequency 
associated  to  this  emission  corresponds  more  or  less  to  the  one  of  the  largest  amplitude  of  the 
eigenmode  at  the  abscissa  from  which  nonlinear  phenomena  start  to  operate.  At  each  time  that 
a  structure  passes  the  throat  (exit  section),  there  is  an  emission  of  a  reflecting  pressure  wave 
(it  has  been  measured  in  the  cold  gas  facilities).  This  wave  arrives  at  the  front  wall  and  there 
is  another  reflecting  wave,  but  now  in  the  downstream  direction.  This  wave  may  reasonably 
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excite  the  intrinsic  instability  mode  (by  the  so-called  receptivity  mechanism)  and  mainly  the  one 
corresponding  to  its  frequency.  This  mode  is  then  fully  favoured  :  it  is  an  intrinsically  amplified 
mode  and  it  starts  with  a  larger  and  larger  amplitude  due  to  the  reflecting  waves.  Rapidly  only 
one  mode  comes  up.  Then,  at  least  in  the  VECLA  and  VALDO  facilities,  there  are  two  cases  : 
either  the  dominant  frequency  corresponds  to  the  one  of  a  cavity  mode  or  not.  The  frequency 
obtained  by  the  instability  and  the  reflecting  waves  strongly  depends  on  the  injection  velocity 
whereas  the  cavity  frequencies  are  obviously  mainly  independent  of  this  injection  velocity.  In 
the  second  case,  the  reflecting  waves  remain  linear  and  the  physics  remains  governed  by  the 
instability  mechanism.  In  the  first  case,  the  amplitudes  are  very  large  and  it  may  be  guessed 
that  the  instability  plays  a  role  only  in  the  transient  which  leads  to  the  excitation  of  the  cavity 
mode. 

•  Finally  when  h  is  large,  the  emission  of  the  large  structure  does  not  occur  any  more  in  the 
booster,  it  moved  downstream,  only  small  fluctuations  (the  intrinsic  instability  modes  them¬ 
selves)  remain  and  do  not  provide  enough  energy  for  creating  reflected  pressure  waves. 

This  scenario  needs  of  course  additional  confirmations  but  it  is  fully  consistent  with  the  experiment, 
the  numerical  simulations  and  the  investigated  nonlinear  stability  approaches. 
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ANNEX  :  Spectral  collocation  method  for  eigenvalue  prob¬ 
lem 


Description  of  the  method 

Several  methods  may  be  used  for  solving  an  eigenvalue  problem  or  more  generally  any  differential 
equations  system.  In  the  case  of  the  plane  Taylor  flow,  expressed  in  terms  of  stream  function, 
equation  (2.9),  the  simplest  method  is  probably  the  spectral  collocation  method.  In  this  annex, 
only  a  short  view  of  it  is  given,  more  details  together  with  other  methods  can  be  found  in  [33], 
an  application  of  the  spectral  collocation  technique  in  a  particular  stability  problem  is  described 
in  [34]. 

The  goal  is  hence  to  solve  the  differential  problem  constituted  by  equation  (2.9)  and  bound¬ 
ary  conditions  (2.10).  The  variable  y  varies  in  [—1,  +1].  Some  definitions  are  necessary,  let  be  first 
Tn  the  nth  Chebychev’s  polynomial  of  degree  n.  Then,  an  integer  N  is  assumed  to  be  chosen,  it 
characterises  the  chosen  refinement  for  solving  the  differential  problem,  as  it  will  be  explained  just 
below.  Let  be  : 


c  I  n3 

il  =  C0S  (  jyT 


j  =  0,...,N 


N  +  1  points,  which  are  called  the  ones  of  “Gauss-Lobatto” .  The  unknown  amplitude  function, 
which  will  be  noted  %f>  in  this  annex,  is  approximated  by  the  polynomial  : 

N 


3=0 


where  A  stands  for 


A  = 


1 zl 


(-i  y+1 


T'n(Z) 

N2a 


with  T'n  the  derivative  of  the  Vth  Chebychev’s  polynomial.  The  N  +  1  discrete  values 

i’j  =  i’iQ 

are  explicitly  the  unknown  of  the  (discretised)  problem.  This  is  the  reason  why  N  describes  the 
used  accuracy  for  computing  the  unknown.  As  ijj  is  a  solution  of  differential  problem,  we  have  now 
to  link  between  the  derivative  ,0/  and  the  function  i/j  itself.  For  this  purpose,  let  be  : 

Cj  (-l)fc+J 


Ejk  — 


£,j  £ k 


Enn  - 


2(1  -  tf) 


3  ±  k 


j  =  1, . . . ,  .ZV  —  1 


E00  =  — Enn  = 


2N2  + 1 
6 


defining  a  square  matrix  E  with  the  coefficients  Cj  given  by 


co  =  cat  =  2 


Cj  =  1 


j  =  1, . . .  ,N  —  1 


Then,  it  can  be  proved  that  the  derivative  exactly  corresponds  for  the  discrete  values  to  the  matrix 
multiplication  by  E  : 


N 


=  ^2Ejki>k 

?  k=0 
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Typical  stability  code 

Using  the  previous  matrix  multiplication  as  derivative  operator  leads  to  discretise  the  linearised 
stability  equations  into  an  algebraic  problem  for  X  the  vector  of  the  unknowns  ipj.  Then,  there 
are  usually  two  approaches,  either  the  complete  spectrum  of  the  discretised  problem  is  looked  for 
or  by  a  shooting  method  a  unique  value  of  the  spectrum  is  looked  for. 

Spectrum.  In  this  case,  if  the  unknown  is  the  complex  number  u>,  the  discretised  problem  is 
written  under  the  form  : 

A.X  =  wB.X 

with  A  and  B  two  matrices  coming  from  derivation  and  intervention  of  the  mean  flow.  All  the 
lines  of  the  matrices  correspond  to  the  discretisation  of  the  continuous  equations  at  the  collocation 
points,  for  example  ordered  from  -1  up  to  +1.  Then  the  four  homogeneous  boundary  conditions  are 
included  by  substituting  the  first  line  by  the  condition  0(— 1)  =  0,  the  second  line  by  0'(— 1)  =  0 
and  the  two  last  lines  in  the  same  way  (by  the  boundary  conditions  in  +1).  Thus  the  previous 
algebraic  problem  is  replaced  by  : 

A.X  =  ujB.X 

The  last  step  consists  in  determining  the  eigenvalues  to  by  an  appropriate  method. 

Shooting  method.  In  this  case,  if  the  unknown  is  the  complex  number  a,  the  discretised  problem 
is  written  under  the  form  : 

A(a).X  =  0 

As  in  the  previous  case,  each  lines  expresses  the  continuous  equation  written  in  the  corresponding 
collocation  point.  The  idea  in  this  case  is  as  before  to  introduce  the  boundary  conditions,  but  only 
three.  In  addition,  a  non  homogenous  one  is  introduced,  for  example  a  condition  on  the  pressure 
(or  the  second  derivative  of  0).  This  non  homogeneous  condition  acts  as  a  normalisation  condition. 
The  previous  problem  is  then  replaced  by  : 

A(a).X  =  b 

The  next  step  consists  in  inverting  A(a)  and  thus  in  determining  X.  Then,  the  omitted  boundary 
condition  is  examined  and  in  fact  all  the  described  procedure  is  included  into  a  loop  :  a  is  searched 
until  the  omitted  boundary  condition  is  satisfied. 

Code  written  in  Matlab 

A  small  program  written  using  the  commercial  software  Matlab  is  given  below.  The  five  routines 
are  organised  as  described  in  figure  4.5.  The  main  program  is  called  “OrrSom”.  The  initialisation 
defines  the  collocation  points  and  the  matrix  used  for  the  derivation.  Then,  it  is  possible  to 
compute  the  whole  spectrum  of  the  linearised  operator,  this  is  achieved  in  the  routine  “spectre”. 
Alternatively,  the  shooting  method  can  be  used,  the  initial  guess  is  given  in  the  main  program, 
the  routine  “balai”  defines  possible  parametric  computations  (by  varying  regularly  one  parameter 
among  the  Reynolds  number,  the  frequency  and  the  abscissa).  This  routine  uses  the  one  named 
“Newton”  whose  goal  is  to  determine  the  eigenvalue  by  a  standard  Newton  iterative  procedure. 
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Figure  4.5:  Representation  of  five  Matlab  routines  for  solving  a  stability  problem  with  the  spectral 
collocation  method  either  with  a  shooting  method  (upper  line  “Balai  -  Newton  -  Resol”)  or  by 
determining  the  whole  spectrum  (lower  line  “Spectre”) 


Each  resolution  of  the  stability  equation  called  by  “Newton”  is  performed  thanks  to  the  routine 
“Resol”. 


Main  program  :  orrsom.m 


Z 

Z  Programme  principal  ORRSOM.m 

Z 

’/,  resolution  du  probleme  de  stability  lineaire 

’/.  pour  l’ecoulement  dans  un  conduit  plan  a  parois  debitantes . 

Z  perturbation  en  forme  de  mode  normal, 

Z  formulation  en  fonction  de  courant 

Z 

Z 

global  npol  xi  e  id  alpha  omega  reynolds  csol  type_calcul  xloc; 

z 

Z  valeurs  de  depart 
Z2.48746e-001  9.70805e-001 
reynolds  =  1000; 
xloc  =  10; 
omega  =  31.6; 

alpha  =  complex(4. 0158,-0. 373) ; 

Z 

type_calcul  =  ’spectre’; 

Z 

initialisation; 

Z 

switch  lower (type_calcul) 
case  ’balai’ 
balai ; 

case  ’spectre’ 
spectre ; 
otherwise 

disp( ’ Orr-Som,  non  implements ’) 

end 
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Subroutine  :  initialisation. m 

•/. 

'/,  Sous-programme  initialisation. m 

•/. 

'/,  Initialisation  dans  la  methode  de  collocation  spectrale 
'/,  npol  est  le  nombre  de  points  de  collocation 

"/,  e  est  la  matrice  de  derivation 

'/,  id  est  la  matrice  identite 

"/,  xi  est  le  vecteur  des  points  de  collocation 

•/. 

global  npol  xi  e  c  id; 

•/. 

'/,  Preliminaires 

•/. 

npol  =  100; 
iaff  =  0; 

xi  =  zeros (npol , 1) ;  c  =  zeros (npol , 1) ;  id  =  zeros (npol ,npol) ;  e  =  zeros (npol, npol) ; 

7. 

"/,  Calcul  du  vecteur  xi  et  du  vecteur  technique  c 

t 

for  j  =  0:npol 

jl  =  j+i; 

xi(jl,l)  =  cos(pi*j/npol)  ;  c(jl,l)  =  1; 

end 

id  =  diag(c , 0) ; 
c (1 , 1)=2 ;  c (npol+1 , 1)=2 ; 

/ 

*/.  Calcul  de  la  matrice  e 

•/. 

for  j  =  0:npol 

ji  =  j+i; 

for  k  =  0:npol 
kl  =  k+1 ; 
if  (k~=j) 

e(  j  1  ,kl)  =  c(jl,l)*(-l)~  (k+j  )  /  (c(kl,l)*(xi(jl,  1)  -xi  (kl ,  1)  )  )  ; 

else 

if  (j~=0  &  j"=npol) 

e ( j 1 ,  j 1)  =  -xi ( j 1 , 1) / (2* (1-xi ( j 1 , 1) *xi ( j 1 , 1) ) ) ; 

end 

end 

end 

end 

e(l,l)  =  (2*npol*npol+l)/6;  e (npol+1 , npol+1)  =  -(2*npol*npol+l)/6; 

Subroutine  :  spectre. m 

7. 
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'/,  Sous-programme  spectre. m 

•/. 

"/,  Determination  du  spectre  de  l’operateur  d’Orr-Sommerfeld 
'/,  on  ecrit  le  probleme  sous  la  forme  : 

"/,  matl.v  =  omega  .  mat2.v 

l 

'/,  Parametres 

/ 


global  npol  xi  e  id  alpha  omega  reynolds  csol  type_calcul  vp  xloc; 

/ 

'/,  initialisation 

•/. 

e2  =  zeros (npol , npol) ;  matl  =  zeros (npol, npol) ;  mat2  =  zeros (npol , npol)  ; 
ci  =  complex(0, 1) ;  e2=e*e  ;  e3=e2*e  ; 

/ 

"/,  ecoulement  de  base 

/ 


u  =  pi*xloc*cos(pi*xi/2)/2;  ddu  =  e2*u; 
v  =  -sin(pi*xi/2) ;  ddv  =  e2*v; 

•/. 


"/,  remplissage  des  matrices  matl  et  mat2 

•/. 


for  i  =  0:npol 
il  =  i+1 ; 
for  j=0:npol 
jl=J+l; 
matl (il , j 1) 
matl (il , j 1) 
matl (il , j 1) 
matl (il , j 1) 

end 


-ci*alpha* (ddu(il)+alpha*alpha*u(il) )*id(il,jl) ; 
matl (il , j l)+ci*alpha*u(il) *e2(il , j 1) ; 
matl (il , j 1)- (alpha*alpha*v(il)+ddv(il) )*e(il,jl); 
matl(il,jl)  +  v(il) *e3(il , j 1) ; 


end 

zl  =-l/reynolds ;  z2  =  2*alpha*alpha/reynolds  ;  z3  =  -alpha. “4/reynolds ; 
matl  =  matl  +  zl*e2*e2  +  z2*e2  +  z3*id  ; 
mat2  =  -ci* (alpha*alpha*id-e2)  ; 

•/. 


'/,  Conditions  aux  limites 


•/. 


for  i  =  0:npol 
il=i+l; 

matl(l,il)  =  0;mat2(l,il)  =  0; 

matl(2,il)  =  e (1 , il) ;mat2(2 , il)  =  0; 

matl (npol , il)  =  e (npol+1 , il) ;  mat2(npol , il)  =  0; 

matl (npol+1 , il)  =  0;  mat2(npol+l , il)  =  0; 

end 

matl (npol+1 , npol+1)  =  1;  matl (1,1)  =  1; 

7. 
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"/,  calcul  des  valeurs  propres  (generalisees) 

l 

bmat  =  eye (size (mat  1) ) ; 

1  =  eig(matl ,mat2) ; 

•/. 

"/,  representation  graphique 

/ 

plot (real (1) , imag(l) , ’ob’ , ’MarkerSize 1 , 6 , ’MarkerFaceColor ’  ,  ’b’ ) 
axis ( [0  100  -10  5]) 

Subroutine  :  balai.m 

7. 

"/,  Sous-programme  balai.m 

7. 

"/,  calcul  de  valeurs  propres  avec  Newton  par  variation  de  l’un  des  parametres 

7. 

'/,  Ce  parametre  peut  etre  : 

"/,  le  nombre  de  Reynolds 

'/,  la  position  en  abscisse 

"/,  la  frequence 

•/. 

7. 

global  alpha  omega  reynolds  xloc; 

7. 

"/.  Valeurs  initiales  pour  la  variation 

•/. 

type_balai  =  ’reynolds’; 
nparam  =  1;  dparam  =  100; 

•/. 

"/,  Lecture  du  fichier  ou  sont  ecrites  les  valeurs  precedentes 

•/. 

fid  =  fopen( ’C:\GREG\ONERA\SPADA\StabMatlab\varirey. txt ’, ’a+’); 

•/. 

[vect,  count]  =  f  scanf  (f  id, ’’/,f  ‘/,f  "/0e  ”/0e  '/,e’,inf); 
nligne  =  count/5;  indice  =  5* (nligne-1) ; 
reynolds  =  vect (indice+1) ;  xloc  =  vect (indice+2) ; 
omega  =  vect (indice+3) ; 

alpha  =  complex(vect (indice+4) , vect (indice+5) ) ; 

•/. 

"/,  Valeurs  pour  demarrer  le  calcul 

•/. 

alphainit  =  alpha;  omegainit  =  omega; 
reynoldsinit  =  reynolds;  xlocinit  =  xloc; 

7. 

"/,  Variation  suivant  le  parametre  choisi 

•/. 

switch  lower (type_balai) 
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case  ’frequence’ 

domega  =  dparam;  dreynolds  =  0;  dxloc  =  0; 
case  ’abscisse’ 

domega  =  0;  dreynolds  =  0;  dxloc  =  dparam; 
case  ’reynolds’ 

domega  =  0;  dreynolds  =  dparam;  dxloc  =  0; 
otherwise 

disp(’non  implements  (ou  erreur)  dans  balai’) 

end 

7. 

’/,  Boucle  pour  la  variation  du  parametre  choisi 

7. 


for  iparam  =  l:nparam 

omega  =  omegainit  +  ( iparam- 1 ) *domega; 
reynolds  =  reynoldsinit  +  ( iparam- 1 ) *dreynolds ; 
xloc  =  xlocinit  +  ( iparam- 1 ) *dxloc ; 

7. 

I  Estimation  des  valeurs  initiates 

7. 

switch  lower (iparam) 
case  1 

alpha  =  alphainit ; 
case  2 


alpha  =  alphal ; 
case  3 

alpha  =  2*alpha2  -  alphal ; 
otherwise 

alpha  =  3* (alpha3-alpha2)  +  alphal; 

end 

7. 

7.  Appel  a  la  convergence,  methode  de  Newton 

7. 

newton; 

7. 

7.  Stockage  des  au  plus  3  dernieres  valeurs  convergees 

7. 


switch  lower (iparam) 
case  1 


alphal  =  alpha ; 
case  2 


alpha2  =  alpha; 
case  3 


alpha3  =  alpha; 
otherwise 

alphal=alpha2 ;  alpha2=alpha3 ;  alpha3=alpha; 

end 

7. 
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7.  Affichage  ecran  des  resultats 

7. 

f  printf  ( ’  iparam  =  7.d,  R  =  ’/,0.1f,  X  =  7.0.  If,  w  =  "/,0.5e,  ar  =  7.0. 5e,  ai  =  7.0 . 5e\n\n\n’ ,  .  .  . 
iparam, reynolds ,xloc , omega, real (alpha) , imag( alpha) ) 

7. 

7.  Ecriture  sur  fichier  des  resultats 

7. 

f printf ( ’ \n’ ) 

fprintf  (fid, ’7.0.  If  7.0.  If  7.0. 5e  7.0. 5e  7.0. 5e  \n  ’,... 
reynolds ,xloc , omega, real (alpha) , imag(alpha) ) ; 

end 

7. 

status  =  fclose(fid); 


Subroutine  :  newton. m 


7. 

7.  Sous-programme  Newton.m 

7. 

7.  Methode  de  Newton  complexe  pour 

7.  resoudre  csol  =  0,  alpha  est  l’inconnue  complexe 

7. 

global  npol  xi  e  id  alpha  omega  reynolds  csol  type_calcul; 

7. 

7.  Initialisation  Newton 

7. 

dval  =  IE-4; 
err  =  1 ; 
it  =  0; 

alphan  =  alpha  ; 

7. 

7.  Methode  de  Newton  (utilisation  des  relations  de  Cauchy) 

7. 

while  (err  >  IE-08  &  it  <  10) 
it  =  it  +  1 ; 
alpha  =  alphan; 
resol;  csolO  =  csol  ; 
alpha  =  alphan  +  dval ; 
resol;  csolm  =  csol; 

dcsol  =  (csolm-csol0)/dval;  cor  =  csolO/dcsol; 
err  =  abs (cor) /abs (alphan) ; 
alphan  =  alphan  -  cor; 

end 

if  (err  >  IE-06) 

disp( ’probleme  Newton’) 
break 

end 
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Subroutine  :  resol. m 


•/. 

'/,  Sous-programme  resol . 

•/. 

"/,  On  ecrit  Orr-Sommerf eld  sous  la  forme  mat.X  =  b 

•/. 

"/,  e2  :  matrice  de  derivee  seconde  e2  =  e*e 

l 


global  npol  xi  e  id  c  alpha  omega  reynolds  csol  xloc ; 

/ 

'/,  initialisation 

/ 

b  =  zeros (npol , 1) ;  e2  =  zeros (npol ,npol) ;  mat  =  zeros (npol , npol) ; 

e3  =  zeros (npol , npol) ; 

ci  =  complex(0, 1) ;  e2=e*e;  e3=e2*e; 

•/. 

"/,  ecoulement  de  base 

/ 

u  =  pi*xloc*cos(pi*xi/2)/2;  ddu  =  e2*u; 
v  =  -sin(pi*xi/2) ;  ddv  =  e2*v; 

/ 


"/,  remplissage  de  la  matrice  mat 

•/. 


for  i  =  0:npol 
il  =  i+1 ; 
for  j=0:npol 
jl=J+l; 
mat (il , j 1) 
mat (il , j 1) 
mat (il , j 1) 
mat (il , j 1) 
mat (il , j 1) 

end 


-ci*alpha* (ddu(il)+alpha*alpha*u(il) )*id(il,jl) ; 

mat (il , j l)+ci*alpha*u(il) *e2(il , j 1) ; 

mat (il , j 1)- (alpha*alpha*v(il)+ddv(il) )*e(il,jl) ; 

mat(il,jl)  +  v(il) *e3(il , j 1) ; 

mat(il,jl)  -ci*omega*e2(il , j 1)  ; 


end 

zl  =-l/reynolds ;  z2  =  2*alpha*alpha/reynolds  ; 
z3  =  -alpha. “4/reynolds  +  ci*omega*alpha*alpha  ; 
mat  =  mat  +  zl*e2*e2  +  z2*e2  +  z3*id  ; 

7. 

'/,  conditions  aux  limites 

"/,  on  ecrit  psi"(l)  =  1,  psi’(l)  =  0,  psi’(-l)  =  0,  psi(-l)  =  0 

7. 

for  i  =  0:npol 
il=i+l ; 

mat(l,il)  =  e2(l,  il) ; 
mat (2 , il)  =  e(l ,  il)  ; 
mat(npol,il)  =  e(npol+l , il) ; 
mat (npol+1 , il)  =  0; 
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b(il)  =  0; 

end 

mat (npol+1 ,npol+l)  =  1;  b(l)  =  1; 

7. 

"/,  inversion 

•/. 

sol_v  =  mat\b  ; 

•/. 

'/,  si  iaff=0,  c’est  Newton  sur  v(l)  ,  sinon  affichage  de  la  fonction  propre 

•/. 

iaff  =  0; 
if  (iaff  ==  0) 

csol  =  sol_v(l) ; 

else 

sol_dv  =  e*sol_v; 

plot(xi,abs(sol_v) ,xi,abs(sol_dv) , ’LineWidth’ ,2) 

set (gca, ’Def aulttextFontName ’ , ’Times  New  Roman’ ) ; set (0, ’Def aulttextFontSize ’ , 14) ; 
xlabel( ’ {\it  y}’ , ’FontName’ , ’Times  New  Roman’ , ’FontSize ’, 16) ; 
ylabel( ’ {\it  u,  v} ’, ’FontName ’, ’Times  New  Roman’ , ’FontSize 16)  ; 
legend (’u’,’v’,0); 

end 
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